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ABSTRACT 

We perform a semi- automated survey for T912 > 2 Lyman Limit systems (LLSs) in quasar spectra 
from the Sloan Digital Sky Survey, Data Release 7. From a starting sample of 2473 quasars with 
z cm = 3.6 — 4.4, we analyze 469 spectra meeting strict seletion criteria for a total redshift path 
Az = 93.8 and identify 192 intervening systems at zlls > 3.3. The incidence of rg 12 > 2 LLSs 
per unit redshift, i r >i{z), is well described by a single-power law at these redshifts: l T >i{z) — 
Clls[(1 + z)/(l + z*)] 7LLS , with z* = 3.7, C LL s = 1.9 ± 0.2, and 7LLS = 5.2 ± 1.5 (68% c.L). These 
values are systematically lower than previous estimates (especially at z < 4) but are consistent with 
recent measurements of the mean free path to ionizing radiation. Extrapolations of this power-law 
to z = are inconsistent with previous estimations of t{z) at z < 1 and suggest a break at z ~ 2, 
similar to that observed for the Lya forest. Our results also indicate that the systems giving rise 
to LLS absorption decrease by w 50% in comoving number density and/or physical size from z = 4 
to 3.3, perhaps due to an enhanced extragalactic ultraviolet background. The observations place 
an integral constraint on the H I frequency distribution /(A/hi, X) and indicate that the power-law 
slope /3 = dln/(A^Hi, X)/d\n Am is likely shallower than /3 = — 1 at Njh f» 10 18 cm~ 2 . Including 
other constraints on f(Nm,X) from the literature, we infer that (3 is steeper than j3 = —1.7 at 
A/hi ~ 10 15 cm~ 2 , implying at least two inflections in f(Nm,X). We also perform a survey for 
proximate LLSs (PLLSs) and find that £plls(z) is systematically lower (w 25%) than intervening 
systems. Finally, we estimate that systematic effects impose an uncertainty of 10 — 20% in the l(z) 
measurements; these effects may limit the precision of all future surveys. 
Subject headings: absorption lines - intergalactic medium - Lyman limit systems - SDSS 



1. INTRODUCTION 

Studies of hydrogen absorption in the lines of sight to- 
wards distant quasars have served to both define, and 
in recent years bring precision to, our cosmological mod- 
els. The low density, highly ionized Lyman-a forest lines 
(a.k.a. the intergalactic medium, IGM), with H I col- 
umn densities A/hi < 10 17 ' 2 cm -2 , have through their ag- 
gregate statistical properties (e.g. their flux power spec- 
trum, mean flux, and column density distributions) con- 
strained cosmological parameters such as the primordial 
power spectrum and the baryonic mass density and as- 
trophysical parameters like the amplitude of the ioniz- 
ing background (e.g. Rauch 1998; Croft et al. 2002; Mc- 
Donald et al. 2005; Tytler et al. 2004; Faucher-Giguere 
et al. 2008b). The high-density, predominantly neutral 
damped Lya systems (DLAs), with A/hi > 10 20 - 3 cm~ 2 , 
trace the gas which forms stars, and likely represent 
the progenitors of modern-day galaxies (e.g. Wolfe et al. 
1995, 2005; Prochaska & Wolfe 2009). 

The majority of Lyman-a forest lines and the DLAs 
have, through analysis of their Lya lines, precisely mea- 
sured Nm values that permit detailed study of their 
physical properties (e.g. metallicity) . For systems with 
intermediate Am values (« 10 18 cm -2 ), however, Lya 
and most of the Lyman series lines lie on the flat portion 
of the curve-of-growth making the Am value difficult to 
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constrain. On the other hand, these systems are optically 
thick to ionizing radiation and impose a readily identi- 
fied signature in a quasar spectrum at the Lyman limit. 
These so-called Lyman limit systems (LLSs), currently 
the least-well studied of H I absorption systems at high 
redshift, are the focus of this manuscript. 

Historically, the LLSs were among the first class of 
quasar absorption line (QAL) systems to be surveyed 
(Tytler 1982). This is because their spectral signature 
is obvious in low-resolution, low S/N spectra. The prin- 
cipal challenge is that the Lyman limit occurs redward 
of the atmospheric cutoff only for systems with redshifts 
z > 2.6. For lower redshifts, one requires spectrometers 
on space-borne ultraviolet satellites. By the mid 1990's, 
samples of several tens of LLSs were generated spanning 
redshifts < z < 4 (Sargent et al. 1989; Lanzetta 1991; 
Storrie-Lombardi et al. 1994; Stengler-Larrea et al. 1995). 
These results were derived from heterogeneous sets of 
quasars discovered from a combination of color-selection, 
radio detection, and slitless spectroscopic surveys. The 
spectra, too, were acquired with a diverse set of instru- 
mentation and therefore varying S /N and spectral resolu- 
tion mitigating differing sensitivity to the precise optical 
depth at the Lyman limit. Although the results were 
not fully consistent with one another, the general pic- 
ture that resulted was a rapidly evolving population of 
absorption systems reasonably described by a (1 + z) 
power-law. 

Cosmologically, the LLSs contribute much if not most 
of the universe's opacity to ionizing radiation. And, until 
recently, the observed incidence of the LLS provided the 
only direct means of estimating the mean free path 
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at any redshift (e.g. Meiksin & Madau 1993; Madau et al. 
1999; Faucher-Giguere et al. 2008a). In a companion 
paper (Prochaska et al. 2009, ; hereafter PWO09), we 
have presented a new technique to measure A 9 ^ 2 that 
circumvents any knowledge of the LLSs. A more precise 
census of the LLSs will serve as a consistency check for 
this A^ f p calculation, but is unlikely to ever again be 
a competitive approach. Instead, the incidence of LLS 
can be used in combination with estimates of A^p to 
assess the Am frequency distribution for gas with Am ~ 
10 16-18 cm~ 2 , a regime that is very difficult to explore by 
studying individual absorption systems. Surveys of the 
LLSs are also likely to place tight constraints on z ~ 3 
cosmological simulations that include radiative transfer. 

Physically, the nature of systems that give rise to a 
LLS remains an open question. The systems with the 
largest TVhi values (i.e. Am > 10 19 cm~ 2 , the so-called 
the super-LLS or SLLS and DLAs) are likely associated 
with the interstellar medium and outer regions of high z 
galaxies. These high A H i systems, however, are only a 
subset of the LLS population. Unfortunately, a proper 
modeling of the LLSs almost certainly requires careful 
modeling of radiative transfer in cosmological simula- 
tions which has thus far been beyond the scope of mod- 
ern computations in cosmological simulations. Indeed, 
the few studies to date have tended to severely underes- 
timate the incidence of LLS (Katz et al. 1996; Gardner 
et al. 2001, but sec Kohler & Gnedin 2007). In recent 
simulations of high z galaxy formation, however, theo- 
rists have placed great attention on 'streams' of cold gas 
that carry fresh material from the IGM to star-forming 
galaxies (Keres et al. 2005; Dekel et al. 2009). These 
cold streams have relatively large hydrogen surface den- 
sities (Ah ~ 10 20 cm" 2 ) and could therefore produce Ly- 
man limit absorption provided the material has a non- 
negligible neutral fraction. Consequently, an accurate 
census of the LLSs with redshift may directly constrain 
the nature and prevalence of cold streams in the young 
universe. 

A final, yet perhaps most important, motivation for 
studying the LLSs is that these systems may dominate 
the census of metals at all epochs. The majority of LLSs 
are metal-bearing, showing metal-line transitions of com- 
mon low and high-ions (e.g. Prochaska 1999; Prochter 
et al. 2009). Because the estimated ionization corrections 
for LLSs with Am « 10 18 cm -2 is large, observations of 
ions in an LLS likely track only a trace amount of the 
metals actually present in the gas. Lyman limit systems 
may show a wider spread in their ionization and metal 
content relative to the IGM or DLA, further emphasizing 
the need for a robust LLS survey. 

In this paper, we survey the homogeneous dataset 
of quasar spectra from the Sloan Digital Sky Survey 
(SDSS), using all 7 public data releases. Our observa- 
tional analysis aims to produce the most precise mea- 
surement of the LLS incidence paying careful attention 
to systematic biases. The wavelength coverage and data 
quality of the SDSS quasar spectra focus the survey at 
z « 3.5. Future work will depend on follow-up obser- 
vations of well-defined quasar samples at other wave- 
lengths. 

The paper is organized as follows. In § 2, we present a 
set of LLS definitions used throughout the manuscript. 



The selection criteria and data quality of quasars from 
the SDSS database are described in § 3. The procedure 
to model the absorbed quasar continuum is presented in 
§ 4 and the search and characterization of LLSs is de- 
tailed in § 5. The criteria used to measure the survey 
path are described in § 6 and an assessment of system- 
atic error and bias from analysis of mock spectra is pro- 
vided in § 7. § 8 describes the principal results and the 
implications for the IGM and cosmology are discussed 
in § 9. Finally, § 10 presents a summary of the main 
findings. Throughout the paper, we adopt a ACDM cos- 
mology with H — 72 /172 km s _1 Mpc -1 , O m = 0.3, and 
fl\ = 0.7 and report proper lengths unless otherwise in- 
dicated. 

2. LYMAN LIMIT SYSTEM DEFINITIONS 

The photon cross-section of a hydrogen atom at ener- 
gies above the Lyman limit may be approximated by: 

<Jli>> V912) « 6.35 x 10~ 18 ( — ) cm 2 , (1) 

\ ^912 / 

with 

^912 = Egi 2 /h = C/A912 , (2) 

and E gl2 = 1 Ryd. Specifically, v> 912 = 3.29 x 10 15 Hz 
and A 9 i 2 = 91 1.7641 A. This implies an optical depth at 
wavelengths A < A912, 

TLL(A - A912) " l0^Tm- 2 (A^) ' (3) 

where Am is the H I column density. For a gas 'cloud' in- 
tersecting a background source with intrinsic flux Fi nt (A), 
the observed flux F Q b s (A) blueward of the Lyman limit is 

-^obs (A < A912) = F int (A)exp [-tll(A)] (4) 

In what follows, we define a 'standard' Lyman limit 
system to be one where the optical depth at A912 is 
T912 > 2, i.e. A H i > 10 17 ' 5 cm -2 . We refer to these sys- 
tems as the T912 > 2 LLS. This corresponds to greater 
than 85% attenuation of an incident ionizing radiation 
field at v = z/912. By this definition, the class of LLS 
includes systems with 10 20 3 cm~ 2 > A H i > 10 19 cm~ 2 
(the so-called super-LLS or sub-DLAs, hereafter referred 
to as SLLS; e.g. O'Meara et al. 2007) and systems with 
Am > 10 20,3 cm~ 2 (the damped Lya systems, DLAs; 
e.g. Wolfe et al. 2005). In a few cases, we will distin- 
guish between these 'strong' LLSs from those with lower 
Am, referring to the latter as T912 < 10 LLS. We also 
note that our T912 > 2 definition for a LLS differs from 
other works which adopted T912 > 1 or T912 > 1.5. These 
are all observationally-driven, not physically-motivated 
definitions. 

Observationally, the absorption of a background source 
by a Tg 12 > 2 LLS is readily apparent, even in low 
S/N spectra. We define absorbers with rgi2 < 2 (i.e. 
Am < 10 17,5 cm~ 2 ) as the partial Lyman limit systems 
(pLLSs). To survey these systems, one requires higher 
quality spectra or an alternate approach to the analysis. 

We define the redshift of an LLS as 
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where Ag^ marks the observed onset of LL absorption. 
In practice, this is is often estimated from strong Lyman 
series lines (e.g. Lya, Ly/3) that accompany the Lyman 
limit opacity. 

We define the sub-set of LLSs that occur within 
3000km s _1 of the emission redshift of the background 
source as proximate LLSs (PLLSs). We separate the 
analysis of these systems from the rest to investigate 
changes in the incidence of optically thick gas near high z 
quasars due to, e.g. the quasar's radiation field and local 
environment. 

Finally, we define the observable £ T >2(z) as the aver- 
age number of T912 > 2 LLS detected per unit redshift 
at a given redshift. In the previous literature, this quan- 
tity is also expressed as n(z), dN/dz, and dn/dz. For 
comparison with previous results in the literature, we 
also consider £ T >i(z), the number of r 9 i 2 > 1 LLSs de- 
tected per unit redshift. We also attempt to separate 
the contributions to £ T >2(z) from SLLSa ^slls( 2 ), DLAs 
^dla(z), and attribute the remainder to the LLSs with 
10 17 - 5 < N m < 10 19 cm 



"LLS 



(Z). 



3. SDSS QUASAR SAMPLE AND SPECTROSCOPY 

One of the primary objectives of the Sloan Digital Sky 
Survey (SDSS) was to discover ~ 100,000 new quasars 
across the northern sky (York et al. 2000). The strat- 
egy of the SDSS team to achieve this ambitious goal was 
a four-fold process: (i) obtain deep, multi-band images 
across a large area of the sky; (ii) select quasar candidates 
by demanding a point-like, point-spread-function and im- 
posing color criteria that separate the candidates from 
the Galactic stellar locus; (iii) obtain follow-up spectra 
for a magnitude-limited sample with a fiber-fed spectro- 
graph. The details of target selection and quasar com- 
pleteness with redshift is described at length in a series of 
SDSS papers (e.g. Richards et al. 2002), but see Worseck 
& Prochaska (2009) for a new and more accurate analy- 
sis; and (iv) automatically identify quasars and estimate 
their redshifts (z em ) through template fitting to the op- 
tical spectroscopy. 

Of these steps, the second has the greatest impact on 
a survey for high z Lyman limit systems. The key issue 
for our survey is whether the presence of an intervening 
LLS biases the targeting of the background quasar for 
follow-up spectroscopy. In effect, a high z LLS severely 
'reddens' the quasar at the bluest optical wavelengths of 
the SDSS imaging. With this effect in mind, the SDSS 
team imposed cuts on the (u — g) color which better sep- 
arated the quasar locus in color pace from the stellar lo- 
cus. The net effect, however, is to bias the spectroscopic 
follow-up against quasar sightlincs without a foreground 
LLS (PWO09). Our analysis indicates an important bias 
for quasars with z om < 3.6. For this reason, we limit the 
statistical analysis to quasars with z cm > 3.6, but we also 
explore the bias by considering the incidence of LLSs to- 
ward quasars with z om = 3.4 — 3.6. 

The quasar spectra analyzed in this paper were taken 
from the Sloan Digital Sky Survey, Data Release 7 
(Abazajian et al. 2009). We retrieved the 'best' ID spec- 
trum for every source flagged as a QSO or HIZ.QSO. This 
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Fig. 1. — Histograms of (top) the emission redshifts z cm for the 
quasars comprising our survey and (bottom) the measurements 
of S/N^ 12 , the signal-to-noise of the absorbed continuum at the 
Lyman limit of each quasar. The sample is restricted to S/Ng- 12 > 
2. 

totaled 102,418 unique spectra 3 . The SDSS survey cm- 
ploys a fiber-fed, dual-camera spectrometer that provides 
continuous wavelength coverage from A « 3800 — 9200A 
at a spectral resolution of FWHM « 150km s _1 . The 
SDSS team employs a custom, data-reduction pipeline 
that performs sky subtraction using empirical measure- 
ments from fibers placed to avoid objects detected in the 
SDSS images. The majority of data suffer from excessive 
sky noise at long wavelengths (A > 8000A) and the in- 
strument throughput and atmospheric absorption limits 
the sensitivity at the shortest wavelengths (A < 4200 A). 

A survey for H I Lyman limit absorption in quasar 
spectra involves two principle steps. First, one must as- 
sess the flux at wavelengths near the Lyman limit in the 
quasar's rest-frame, A < Agi2(l + z enl ). We discuss our 
procedure for this step in the following section. The sec- 
ond step is to estimate the flux at wavelengths blueward 
of LLS candidates. There are several characteristics of 
the SDSS spectroscopy which negatively affect this es- 
timate. A generic concern is the poorer instrument re- 
sponse at the bluest wavelengths. At the bluest wave- 
lengths, many of the spectra exhibit a very low signal-to- 
noise ratio and yield flux estimates consistent with zero, 
even without an intervening LLS. Therefore, we have lim- 
ited our survey to redshifts zlls > 3.3 corresponding to 

3 The SDSS spectra were processed through our automated al- 
gorithms for finding absorption-line features and damped Lya can- 
didates (Prochaska et al. 2005; Herbert- Fort et al. 2006). Approx- 
imately 10 of the spectra failed to be processed (primarily because 
the SDSS-reported emission redshift is erroneous) and were re- 
moved from any subsequent analysis. 
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A > 3920A. In practice, we restrict the quasar sample to 
objects with z em > 3.6 and perform a search for LLSs at 
all redshifts, but then only analyze absorption systems 
with zlls > 3.3. We also limit the survey to quasars 
with z cm < 5 because the SDSS spectra of higher red- 
shift objects are generally too low S/N to permit a robust 
analysis. Figure la presents histograms of the emission 
redshifts for the statistical sample and the signal-to-noise 
of the absorbed continuum at the Lyman limit (S/Ng 12 ! 
see below for the definition). 

Another difficulty with the SDSS spectra at blue wave- 
lengths is that the two-dimensional spectra of faint 
sources may be improperly traced. On occasion, the ID 
extractions include flux from a neighboring object and 
yield a systematic overestimate of the flux. This effect 
reduces the estimated opacity for a LLS. Another issue, 
especially with a fiber- fed spectrometer, is that the sky 
model is estimated from nearby fibers that are intention- 
ally placed on 'object-free' regions of the sky. Although 
the SDSS project has worked carefully to mitigate the ef- 
fects of variable fiber efficiency, significant misestimates 
of the sky are known to occur. We have identified tens of 
objects where the extracted flux is significantly negative, 
indicating an overestimate in the sky model. This may 
convert a partial LLS (with Tgi2 < 1) into a rgi2 > 2 
LLS. By a similar token, a proper T912 > 2 LLS may ap- 
pear as a pLLS if the sky is underestimated. We proceed 
under the expectation that this effect is nearly random, 
i.e. for every underestimate of the sky there is a corre- 
sponding overestimate, but this has not been rigorously 
established. The SDSS fibers are sufficiently wide (diam- 
eter of 3") that they will occasionally include flux from 
a projected neighbor. These coincident objects may be 
much fainter than the quasar at redder wavelengths, but 
they could contribute all of the flux blueward of a strong 
LLS and lead to an underestimate of the LL opacity 4 . 

4. THE ABSORBED QUASAR CONTINUUM 

Absent any other sources of opacity, one can trivially 
estimate rgi2 from the quasar spectrum through mea- 
surements of the flux both redward and blueward of the 
observed Lyman limit (Equation 4). In practice, how- 
ever, the quasar flux is also attenuated by line opacity 
from the so-called Lya forest (a.k.a., the intergalactic 
medium; IGM). For example, consider a Lyman limit sys- 
tem at zlls = 3.5 intervening a z cm — 4 quasar. The LLS 
attenuates the quasar flux blueward of Agf 2 s = 4103A. 
At this wavelength, the quasar spectrum recorded on 
Earth will also include opacity from the Lya forest at 
= (1 + z L Ls)(A9i2/ALya) - 1, Ly/3 absorption from 
the IGM at z Ly/3 = (1 + ZLLs)(A 9 i2/A L y/3) - 1, etc. It 
is necessary, therefore, to account for these additional 
sources of opacity when estimating Tgi2. 

We can express the observed (rest-frame) quasar flux 
-Fobs in terms of the intrinsic flux (just) redward of the 
Lyman limit Fi nt as 

Fobs (A > A912) = F lnt exp[-T IGM (A)] , (6) 

4 An amusing (and plausible) systematic effect related to this 
is contamination by the light reflected from terrestrial satellites 
crossing the night sky. Even a brief 'exposure' through the 3" 
fiber could dominate the flux at the bluest wavelengths, although 
this should generally be mitigated by the fact that the SDSS team 
acquires 3 unique exposures per target. 
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Fig. 2. — Mean observed quasar spectra from the SDSS- 
DR3 in redshift intervals from top to bottom of z = 
[3.4,3.7); [3.7,4.0); [4.0,4.4) and [4.4,5.0). The spectra have been 
offset by a factor of 1.5 for presentation purposes and the normal- 
ization is not physically meaningful. The increasing impact of the 
Lya forest is readily apparent at wavelengths A ros t < 1200A. At 
Arcst < 920A, we have set the template to a constant value. 

where tigm is the effective opacity of the IGM from Ly- 
man series line-opacity 5 . The LLS introduces an addi- 
tional, continuous opacity blueward of the Lyman limit: 

Fobs (A < A912) = F int exp[-Ti GM (A) - r LL (A)] . (7) 

A precise estimate of T912 = TLL(Agi 2 ), therefore, requires 
an estimation of the absorbed quasar flux not its intrin- 
sic flux. Conveniently, this quantity is the observed flux 
recorded in the spectrum at A > A912. There are still sig- 
nificant challenges because the IGM opacity is stochas- 
tic on both small (individual Lyman lines) and modest 
scales (many 10A) and the intrinsic quasar spectrum 
(both shape and normalization) varies from source to 
source. We now describe an automated procedure used 
to estimate the absorbed continuum from each quasar 
spectrum. 

The traditional method of estimating the quasar con- 
tinuum is to first identify regions of unabsorbed quasar 
flux and then to interpolate a continuum level between 
these regions. For quasars at high redshift, this method 
is particularly error prone, because at wavelengths below 
Lya emission we expect few (and at very high redshift, 
none) of the pixels to be free of absorption from the IGM. 
Moreover, this traditional method frequently requires by- 
hand modification, which is time-intensive and subjective 
to individual biases. Methods do exist to automatically 
generate a quasar continuum from emission-line charac- 
teristics (e.g. Suzuki 2006), but these are designed to in- 
fer the intrinsic quasar spectrum not the IGM-absorbed 
continuum. 

Our approach is to match a template model of the aver- 
age absorbed continuum to each spectrum, allowing for a 
large-scale tilt (i.e. a unique underlying power-law slope) 
and arbitrary normalization. We emphasize that this ap- 
proach is not intended to recover the intrinsic spectral 
energy distribution of the quasars. Our scientific inter- 

5 In the following, we do not explicitly derive the opacity from 
metals in the IGM, but these may be considered included in tiqm- 
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TABLE 1 
SDSS-DR7 QUASAR SAMPLE 



Plate 


MJD 


FiberlD 


Object Name 


2em 


fig a 


A b 


B b 


Scale c 


S/NA 12 


650 


52143 


111 


J000238. 41-101149.8 


3.938 





2.60 


-0.29 


1.11 


0.9 


750 


52235 


608 


J000300. 34+160027.7 


3.698 





16.35 


-3.54 


1.00 


1.5 


650 


52143 


48 


J000303. 34-105150.6 


3.646 





46.77 


-11.08 


0.96 


2.0 


750 


52235 


36 


J000335.21+144743.6 


3.484 





13.96 


-3.23 


1.12 


0.8 


751 


52251 


207 


J000536.38+135949.4 


3.686 





9.96 


-2.00 


1.00 


1.1 


751 


52251 


562 


J000730.82+160732.5 


3.501 





7.35 


-1.26 


1.15 


1.2 


651 


52141 


534 


J001001. 02-090519.1 


3.720 


2 










751 


52251 


39 


J001115.23+144601.8 


4.967 





80.36 


-18.34 


1.00 


3.0 


752 


52251 


378 


J001134.52+155137.3 


4.325 





15.67 


-3.60 


1.00 


1.1 


752 


52251 


204 


J001328. 21+135828.0 


3.575 





16.28 


-2.75 


1.07 


2.5 


752 


52251 


5 


J001747.90+141015.7 


3.955 





21.69 


-5.06 


0.92 


1.0 


753 


52233 


310 


J001813.88+142455.6 


4.235 





22.74 


-5.04 


1.00 


1.7 


753 


52233 


291 


J001820. 71+141851. 5 


3.936 





70.03 


-16.86 


1.32 


3.6 


753 


52233 


391 


J001918.43+150611.3 


4.134 





1.91 


0.13 


1.13 


1.2 


390 


51900 


271 


J001950.05-004040.7 


4.327 





21.55 


-4.96 


0.89 


1.9 



Note. — [The complete version of this table is in the electronic edition of the Journal. The printed 
edition contains only a sample.] 

a Flag indicating: 0=Normal; l=Not at SDSS reported redshift and/or not a quasar (excluded); 2=Too 
low S/N for evaluation (excluded); 3=Strong BAL (excluded) 

b Absorbed continuum fitting parameters of the form: C = A + B log 10 (A/A). We caution the reader 
that these models arc not meant to describe the intrinsic spectral energy distributions of the quasars (see 
the text). 

c Additinal scaling factor imposed by the authors on the best-fit absorbed quasar continuum. 
d Estimate of the signal-to-noise for the absorbed quasar continuum at A = A9i2(l + z C m). 



est, in this manuscript, is to model the absorbed contin- 
uum of a quasar near its Lyman limit with the fewest 
number of parameters. Indeed, for some quasars we de- 
rive power-law slopes that are likely unphysical due to 
the stochastic nature of IGM absorption, the non-linear 
effects of emission lines, intrinsic deviations from a sin- 
gle power-law, and the narrow spectral range considered. 
We caution against using any of the following results for 
studies on the physical properties of high z quasars. 

Our first step is to derive the templates for the average 
absorbed quasar continuum. Because the line-density of 
the IGM and therefore tiqm increase with redshift, we 
perform this analysis in small redshift intervals (Sz = 0.3 
to 0.6; Figure 2). For every quasar within the redshift 
interval, we shift the spectrum to the quasar rest-frame 
using the SDSS reported z em value. Next, we "dctilt" the 
spectrum by removing a power-law shape. This power- 
law is determined as follows. 

We have constructed from the SDSS-DR3 dataset 
an average template spectrum, archived in XIDL 6 as 
"fulLSDSS-LLS.fits" . For each individual quasar spec- 
trum, we sample quasar pixels with wavelengths greater 
than Lya emission, divide by the template spectrum, 
and measure the slope of the resulting spectrum. Af- 
ter aligning the spectra in the quasar rest-frame (near- 
est pixel), we median-combine the data in each redshift 
interval 7 . The resultant template represents the me- 
dian intrinsic quasar continuum modulated by the me- 
dian flux decrement of the IGM. Due to the presence of 
Lyman limit systems, the quasar templates will not be 
useful at wavelengths near and below the rest-frame Ly- 
man limit; this portion of the spectrum, however, can 
be used to constrain the mean free path to ionizing ra- 
diation (PWO09). For wavelengths blueward of 920A, 

6 http:/ /www. ucolick.org/~xavicr/IDL 

7 This stack is not optimal for deriving Ly« forest statistics 
(Dall'Aglio ct al. 2008). 



therefore, we set the template to have the value recorded 
at 920A. Conveniently, there are no strong emission fea- 
tures in the quasar SED at these wavelengths (e.g. Tclfer 
et al. 2002). The resultant template spectra are shown 
in Figure 2. 

With the templates constructed, a model of the ab- 
sorbed continuum for each quasar is determined as fol- 
lows. First, we shift the observed quasar spectrum to 
the rest-frame and divide by the appropriate template 
spectrum (i.e. according to z om ). Second, we sample the 
quasar in the wavelength range 950A < A < 1800A and 
fit a power-law (p(A) = A + B log[A/A]) to the observed 
flux, weighting by the inverse variance array. The emis- 
sion lines in this spectral range may bias the fit, but we 
do review and modify these fits (see below). The prod- 
uct of this power-law with the template, when shifted to 
the observed frame, provides our model for the absorbed 
quasar continuum. The power-law parameters derived in 
this fashion are listed in Table 1. Sample fits are shown 
in Figure 3. With these models, we can calculate the 
ratio of the absorbed continuum to the lcr-error array 
each quasar's Lyman limit, which we denote as S/N^ 12 . 
Figure lb shows the distribution of S/Ng 12 values for the 
statistical survey. 

5. LYMAN LIMIT SYSTEM SEARCH AND 
CHARACTERIZATION 

In the following, we parameterize a LLS by two quanti- 
ties: (1) its absorption redshift zlls and (2) the total H I 
column density Ahi- Although the H I Lyman series lines 
are sensitive to the component structure and the Doppler 
parameters (also known as b- values) of the 'clouds' com- 
prising an LLS, the opacity blueward of A912 is insensi- 
tive to these details. Furthermore, the SDSS spectra are 
generally of too poor quality to constrain such structure 
using the observed Lyman series lines. Therefore, our 
model of an LLS assumes a single cloud with a Doppler 
parameter of b = 30 km s _1 . An implication of this pa- 
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rameterization is that two systems with small redshift 
separation are modeled as a single system with the total 
of the TVhi values. Our tests with mock spectra (§ 7) 
indicate that two absorbers with Sz < 0.1 are often in- 
distinguishable from a single LLS. The survey presented 
here, therefore, refers to LL absorption smoothed over a 
redshift interval of Sz « 0.1. We return to this point in 
our presentation and discussion of the survey results. 

We have developed an algorithm (sdss_findlls) to au- 
tomatically search for and characterize LLS absorption 
in quasar spectra. In brief, the code generates a set of 
model spectra for the line and continuum opacity of a sin- 
gle LLS with redshifts covering z = z — (A0/A912 — 1) 
to (2; om +0.2) where Ao is the starting wavelength of the 
SDSS spectrum and z em is the quasar emission redshift 
reported in the DR7 8 . The grid of models assumes iVui 
column densities logAni = 16.0, 16.2, 16.4, 19.8 and 
a Dopplcr parameter b = 30km s _1 . We implement a 
grid with 0.2 dex spacing in TVhi because very few of the 
spectra have sufficient S/N to provide a more precise es- 
timate. Furthermore, we estimate systematic uncertain- 
ties (e.g. related to continuum placement, sky subtrac- 
tion) to be of this order. These models arc convolved 
with the SDSS instrumental resolution and then applied 
to the absorbed quasar continuum. Finally, the code 
constructs a x 2 grid in ^lls an d Am space, identifies 
the minimum x 2 , and records the 'best-fit' values. This 
approach differs from previous methods which focused 
solely on the Lyman limit (e.g. Storrie-Lombardi et al. 
1994) or relied on 'by-eye' analysis (Lanzetta 1991). In 
general, the spectra provide very little constraint on Ahi 
for values exceeding 10 17 5 cm~ 2 until the Lya profile be- 
comes damped (e.g. Prochaska et al. 2005). Therefore, 
we report lower limits to Am for any LLS with 7912 > 2. 

For sightlincs with a single LLS having Am > 
10 17 2 cm -2 and good S/Ng\ 2 (i.e. greater than 5 pix -1 
at A|™ 2 = A912 x [1 + Zem]), w e find that the automated 
algorithm is highly successful on its own. In practice, 
however, there are several aspects of the data and anal- 
ysis that require visual inspection of the spectra and in- 
teractive modification to the model: First, many spec- 
tra have such low S/N that zlls and Am are poorly 
determined. In these local minimum in \ 2 can 

occur which gives a mis-estimate of these quantities. Sec- 
ond, the presence of multiple absorbers along the sight- 
line (e.g. one or more pLLSs with a lower redshift LLS) 
gives a spectrum that cannot be well modeled by a single 
LLS. Third, a non-negligible number of the spectra re- 
trieved from the SDSS database purported to be high z 
quasars are either at a lower redshift or are another class 
of astronomical object altogether. Fourth, we found that 
half of the absorbed continuum models required scaling 
to higher or lower value by greater than 10%. Finally, 
we prefer to avoid quasars with with strong broad ab- 
sorption line (BAL) or associated systems to focus the 
analysis on intervening LLSs. 

Given the above complications to an automated anal- 
ysis, we built a graphical user interface (GUI) within 
the IDL software package (sdss-chklls; bundled within 
XIDL) that inputs the data and best-fit LLS model 
for each object. Two of the authors (JXP and JMO) 

8 Modified by the analysis of Prochaska et al. (2008) where ap- 
plicable for quasars from the DR5 release. 



TABLE 2 

SDSS-DR7 INTERVENING r 912 > 2 LLS SURVEY 



Quasar z cm z s / art z st ' art z LLS 



z cm < 3.6 



J001328. 21+135828.0 


3.575 


3.443 








JUIO 1 . do — uiuozy.D 


O. OD^t 


*Z ^87 
0.00 1 


*Z ^S7 
0.00 1 


O 
O 


^87 


J073947.17+445236.7 


3.575 


3.300 








J074914.13+305605.8 


3.436 


3.300 


3.300 






J075303. 34+423130.7 


3.590 


3.300 


3.300 






J075859. 81+165811. 8 


3.439 


3.366 




3 


366 


J080025. 10+441723.1 


3.554 


3.300 


3.362 






J080525. 53+123438.7 


3.425 


3.300 


3.300 










..6 








J001115. 23+144601. 8 


4.967 


4.567 








J001820.71+141851.5 


3.936 


3.536 


3.596 






J004219. 74-102009.4 


3.880 


3.633 


3.633 


3 


633 


J010619.24+004823.3 


4.449 


4.049 


4.049 






J011351.96-093551.0 


3.668 


3.615 


3.615 


3 


615 


J012403. 77+004432.7 


3.834 


3.434 


3.434 






J015048. 82+004126.2 


3.702 


3.302 


3.302 






J015339. 61-001104.8 


4.194 


3.879 




3 


879 



Note. — [The complete version of this table is in the 
electronic edition of the Journal. The printed edition con- 
tains only a sample.] 



used this GUI to validate and/or modify all of the mod- 
els. These authors flagged erroneous spectra (159 ex- 
amples), strong BAL or associated absorption (quasars 
showing very strong C IV, N V, and O VI absorption; 290 
quasars), or data with such low S/N that any analysis 
was deemed impossible (114 spectra). For the remainder 
of sightlincs, the authors could modify the continuum 
(via a multiplicative scalar; Table 1) and/or change the 
model of LLS absorption (i.e. zlls, Am)- This includes 
absorption due to candidate pLLSs. In many cases, zlls 
was modified to correspond to the strongest, local Lya 
absorption line at A = (1 + zlls) x 1215.67A, especially 
for those systems that also showed absorption at the ex- 
pected wavelength for Ly/3. 

After every sightline was analyzed in this manner, the 
results from the two authors were compared to assess 
consistency. Roughly half of the spectra were reviewed 
because of conflicts in the models. The majority of these 
were associated with the absorbed continuum placement 
(typically offsets of 5 — 10%) which implied differences in 
the search path of \Az\ > 0.1 (see § 6). In the majority 
of these cases, we simply averaged the two estimations 
of the continuum. The second most frequent conflict 
was on the definition of strong BAL absorption, primar- 
ily because we did not adopt uniform or strict criteria. 
In most cases, we conservatively excluded the sightline. 
There were also « 100 cases where one author estimated 
log -/Vhi = 17.4 when the other estimated log Am = 17.6, 
i.e., straddling the Tgi2 = 2 boundary that defines our 
LLS search. These were especially scrutinized for the 
presence of higher-order Lyman series lines. Where nec- 
essary, the final best-estimate for Am was deferred to 
the third author (GW). 

Tables 2 and 3 list the set of Lyman limit systems iden- 
tified in the SDSS-DR7 for systems that (respectively) 
influence our statistical analysis and otherwise. For each 
system, we list our best estimate for zlls and the Am 
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TABLE 3 

SDSS-DR7 NON-STATISTICAL LLS SYSTEMS AND 
CANDIDATES 



Quasar 






ZLLS 


loe Nvn a 


S/N 6 


J000238.41-101149 


8 


3 


938 


3 


809 


17.2 


0.7 


J000300.34+160027 


7 


3 


698 


3 


570 


17.2 


1.3 


J000303.34-105150 


6 


3 


646 


3 


467 


17.6 


1.5 


J000335.21+144743 


6 


3 


484 


3 


498 


16.6 


0.9 


J000536.38+135949 


4 


3 


686 


3 


580 


19.0 


1.1 


J000730.82+160732 


5 


3 


501 


3 


511 


19.8 


1.3 


J001115.23+144601 


8 


4 


967 


3 


995 


17.8 


2.6 


J001134.52+155137 


3 


4 


325 


1 


348 


19.8 


1.1 


J001328.21+135828 





3 


575 


3 


282 


19.8 


1.7 


J001747.90+141015 


7 


3 


955 


3 


925 


17.4 


1.0 


J001813.88+142455 


6 


4 


235 


i 


151 


17.2 


1.6 


J001820.71+141851 


5 


3 


936 


3 


456 


17.2 


2.4 


J001918.43+150611 


3 


4 


134 


4 


053 


17.8 


1.2 


J001950.05-004040 


7 


4 


327 


4 


047 


18.4 


1.8 


J002120.05+155125 


7 


3 


698 


3 


671 


17.8 


2.0 


J002614.69+143105 


2 


3 


973 


3 


895 


19.8 


1.3 



Note. — [The complete version of this table is in the 
electronic edition of the Journal. The printed edition con- 
tains only a sample.] 

a The Nhi values listed serve as a rough estimate. Typical 
uncertainties for systems with Ahi < 10 17 ' 5 cm -2 are at 
least 0.2 dex. 

b Estimatc of the S/N of the absorbed continuum (ignoring 
the effects of pLLS) at A = Agi2(l + zlls)- 

value. Because the data offer minimal constraint on Ahi 
for values exceeding 10 175 cm -2 , in the statistical sample 
we only report a lower limit to Art for all systems with 
best-estimates of Ahi > 10 175 cm~ 2 . For the pLLSs, 
the typical la uncertainty is approximately 0.2 dex for 
Ahi > 10 17 cm~ 2 and is dominated by systematic error 
in the continuum placement and the stochastic nature of 
IGM opacity. As noted above, zlls was frequently mod- 
ified in the interactive analysis to correspond to strong 
Lya and Ly/3 lines. Our analysis of mock spectra indi- 
cate typical redshift uncertainties of a(z) < 0.02 with 
occasional 'catastrophic failures' due to line-blending or 
spurious spectra (§ 7). No attempt was made to improve 
this estimate by searching for the presence of metal-line 
absorption (e.g. C IV) outside the Lya forest. There is a 
tendency, both in the automated algorithm and in inter- 
active analysis, to set zlls a t the highest value possible 
that can be accommodated by the data. To this extent, 
we suspect that there is a modest bias in our zlls values 
to slightly higher redshifts (less than 0.01 on average, but 
with the occasional large offset). 

Figure 3 shows a representative sample of four quasar 
spectra, zoomed into the region blueward of Ly/3, with 
the absorbed continuum and LLS models indicated. We 
provide snapshots of the LL region for all quasar spectra 
in the statistical sample online 9 . 



6. SURVEY PATH 

Analogous to galaxy surveys where one defines a search 
volume based on the depth of imaging and spectroscopic 
follow-up, measurements of the incidence of quasar ab- 
sorption line systems requires an estimate of the total 
spectral path sensitive to a robust search. This is gener- 
ally referred to as the redshift path covered (by translat- 
ing observed wavelength into redshift). For the survey 

9 http:/ /www. ucolick.org/~xavier/SDSSLLS 




4500 5000 5500 6000 

Wavelength (Ang> 

Fig. 3. — Representative SDSS-DR7 spectra of quasars from our 
statistical survey. Each example shows the absorbed continuum 
model (dark blue) and the the model that includes LL absorption 
(cyan). The smoothed ler error array is given by the orange curve. 
For this presentation, we have not shown the Lyman series in the 
model. The vertical dotted line indicates the Lyman limit wave- 
length in the quasar rest-frame, i.e. at A^j = Agi2(l + z om ). The 
red bins indicate the mean flux in intervals of 100 pixels. The top 
panel shows a typical example of a T912 > 2 LLS in data that just 
satisfies our S/N^ 12 criterion. The next lower panel shows an ex- 
ample where we identify a pLLS candidate and then a T912 > 2 
LLS. The next lower panel shows a high S/N^ 12 example whose 
quasar spectrum rapidly declines with decreasing wavelength be- 
low Ag" 2 - This absorbed continuum was modelled by including a 
pLLS candidate. The lowest panel shows a low S/N^ 12 example 
where a T912 > 2 LLS is apparent. 

of LLSs, we have adopted the following criteria for in- 
cluding spectral regions in the search. These are based 
on our automated and interactive analysis of the SDSS 
spectra, our modeling of Keck/LRIS follow-up spectra, 
and our analysis of simulated spectra (§ 7, Appendix A): 

1. The search path will begin at a minimum redshift 

Of Z s tart > 3.3 

2. For the intervening LLS sample, the search 
path ends at the redshift z cn( \ corresponding to 
3000 km s _1 (relativistic) blueward of the quasar 
redshift z cm . 

3. The absorbed continuum flux must exceed twice 
the estimated error array, i.e. S/N A 12 . 

4. The search path begins at a maximum offset of 
Sz = z cm - z start < 0.4. 

The first criterion is motivated by the starting wave- 
length of the SDSS spectra (A ~ 3800A) and the poorer 
quality of the data at the bluest wavelengths. For red- 
shifts less than 3.3, there is insufficient spectral coverage 
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and/or data quality to confidently assess the presence of 
an LLS. The second criterion minimizes the influence of 
the quasar and its environment on the analysis. This 
criterion is relaxed in the study of PLLSs. 

The third criterion is the most subjective, yet impor- 
tant, for setting the redshift survey path of each quasar. 
Algorithmically, we impose this constraint by identifying 
the first pixel blueward of Ag 1 -^ = ^912(1 + z cm ) where 
our model of the absorbed quasar continuum falls below 
twice the median-smoothed (15 pixels) la error array. 

This pixel defines the starting redshift z^/^ t ~ 2 corre- 
sponding to a S/Ng 12 =2 limit. If the first pixel blueward 
of Ag™ 2 does not satisfy the S/N threshold, the quasar 
has zero redshift path, i.e. z s tart = ^end- One can, of 
course, define starting redshifts corresponding to higher 
(or lower) S/Ng\ 2 limits; indeed, our fiducial choice of 
S/Ng 12 =2 should be considered arbitrary, although it is 
guided by our analysis of real and simulated spectra. 
And to avoid a systematic bias associated with pLLSs 
(see § 7.3), one must choose the S/Ng 12 criterion to be 
sufficiently high to discover r gl2 > 2 LLS even in the 
presence of a pLLS 10 . We investigate the impact of this 
choice on our results later in the manuscript. Finally, our 
fourth criterion is imposed to mitigate the cumulative ef- 
fects that pLLSs can have on our ability to detect LLS 
with Tgi 2 > 2. That is, trials with mock spectra (§ 7) in- 
dicate that mulptiple pLLSs along a single sightline may 
prevent the detection of a rgi 2 > 2 LLS and that this 
bias is minizmized provided 5z < 0.4. Furthermore, we 
find that the extrapolation of the absorbed continuum 
from the quasar's Lyman limit often is a poor model for 
5z > 0.4. 

The starting redshift is further modified by the pres- 
ence of Lyman limit absorption. In the case of a rgi 2 > 2 
absorber, the quasar flux is severely depressed below 
A = Agi 2 (l + zlls) and we terminate the search path at 
this wavelength. Specifically, this implies 2: sta rt > ^lls 
for all sightlines with an rgi 2 > 2 LLS. For sightlines 
where one or more pLLSs are identified 11 , we had origi- 
nally intended to terminate the search once the absorbed 
quasar continuum convolved with the pLLS absorption 
failed to satisfy the S/N threshold. In our analysis of 
mock spectra, however, we found that this introduces 
a "pLLS-bias" where the incidence of Tgi 2 > 2 LLS is 
overestimated (§ 7.3). In part, our S/Ng\ 2 =2 threshold 
is chosen so that one can robustly search for rgi 2 > 2 
LLSs even along sightlines where one or more pLLSs are 
present. 

Table 2 presents the list of quasars in SDSS-DR7 that 
(i) have 3.6 < z cm < 5, (ii) were not identified to exhibit 

strong BAL signatures, and (iii) have 3.3 < z^/^ t ~ 2 < 
z cm . There are 469 quasars satisfying these criteria. For 
each sightline, we list the starting redshifts for S/Ng\ 2 =2 
and 3 limits. We report results for these two values of the 
S/N threshold to search for data-quality biases. Table 2 

S/N A — 2 

also lists all Lyman limit systems with zlls > ^start 

10 Contrary to some of the earliest work on LLSs (e.g. Tytler 
1982), we do not use survival statistics to estimate the number of 
LLSs at redshifts below any observed LLS. 

11 Again, all pLLSs should be considered candidates. Some of 
the spectral features modeled as pLLSs are instead due to unusual 
variations in the absorbed quasar continuum (e.g. highly reddened 
quasars). 
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Fig. 4. — Summary of the redshift path surveyed in the SDSS- 
DR7 for T912 > 2 LLS absorption assuming S/N^ 12 thresholds of 
2 and 3 (solid curves). Specifically, g(z) represents the number of 
unique quasars in the SDSS that provide a search for LLSs over 
the interval dz = 0.001 at redshift 2. The dotted curve, meanwhile, 
represents the same quantity but for proximate LLSs (PLLSs; sys- 
tems within 8v = 3000km s — 1 of the quasar) and corresponds to 
a S/N* 12 =2 threshold. This curves fall rapidly for z < 3.6 be- 
cause we limit the analysis to quasars with z om > 3.6. The non- 
monotonic nature of g(z) is due primarily to the presence of LLSs 
along the quasar sightlines which stochastically truncate the search 
path. 

We do list these quantities for quasars with z cm < 3.6, 
but these were not included in our final statistical analy- 
sis because of the bias related to SDSS targeting criteria 
previously mentioned by PWO09. 

Using the values presented in Table 2, it is straightfor- 
ward to calculate z cn d and the redshift search path for 
each quasar: Azi = z cn d — z stai t- For a given S/Ng\ 2 
limit, the total search path for the full dataset is 

AZTOT = AZ * ■ ( 8 ) 

We calculate Az^q^~ 2 = 96. To maintain a homoge- 
neous sample and set of search criteria, we do not include 
previous studies in our analysis. We compare against 
previous results in § 9. 

In Figure 4, we present the sensitivity function g(z) 
which expresses the number of SDSS quasars at redshift 
z where a robust search for LLSs is possible. The sev- 
eral solid curves represent differing S /Ng\ 2 limits for the 
survey. We also present the sensitivity function for prox- 
imate LLSs (PLLSs) which corresponds to the sample of 
quasar spectra that satisfy the S/N gl2 > 2 criterion at 
3000 km s _1 blueward of z cm . 

With the definition of the search path and the identi- 
fication of the Lyman limit systems along each sightline 
(§ 5), it is straightforward to calculate the incidence of 
intervening r g i 2 > 2 LLS per redshift interval, ^ r > 2 (z). 
The standard estimator is to compare the total number 
of LLSs against the total survey path in discrete redshift 
intervals. We will return to evaluate i r >i{z) and discuss 
the values after exploring several sources of systematic 
error. 

7. MOCK SPECTRA AND SYSTEMATIC ERRORS 
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With a survey the size of SDSS, one can quickly reduce 
the statistical noise in measurements to very small levels. 
In this regime, one must carefully assess all sources of 
systematic error as these may dominate the measurement 
uncertainty. To this end, we have conducted a range of 
tests with mock spectra as described in this section. The 
casual reader may wish to skip to the summary of this 
section (§ 7.6) 

7.1. Mock Spectra Construction 

We generated a set of 800 mock SDSS quasar spectra 
and analyzed them in the same way as the real data in 
order to assess bias and completeness in our LLS sur- 
vey. The H I forest absorption spectra were generated 
via a Monte Carlo routine similar to the one described 
in Dall'Aglio et al. (2008), assuming that the Lya for- 
est is well characterized by three independent distribu- 
tions: (i) the Lya line-incidence £ a (z) oc (1 + z) 1 , (ii) 
the H I column density distribution f(Nm) oc iV^j, and 
(iii) the Dopplcr parameter distribution parameterized 
as f(b) oc fo- 5 cxp[-^/&4] (Hui & Rutledge 1999). Each 
simulated line of sight was filled with H I Lya absorp- 
tion lines at 2 < z < 4.6 until the H I Lya effective 
optical depth was consistent with Fauchcr-Gigucre ct al. 
(2008b), both in normalization and redshift evolution. 
We did not model the z ~ 3.2 dip in the effective op- 
tical depth measured by Faucher-Giguere et al., and in- 
stead adopted a simple power-law = 0.0011(l+z) 4 23 . 
If the number of lines in a given patch of the forest is 
Poisson-distributed, a power-law line density evolution 
£ a (z) oc (1+z) 7 yields a power-law effective optical depth 
evolution t" s oc (1 + z) 7+1 (Zuo 1993). The column den- 
sity distribution was modeled with a single power-law 
index [3 = -1.5 for 12 < log(iV HI ) < 19, but with a 
0.5 dex break at log(A^ni) = 14.5 in order to account for 
the dearth of high column density lines, consistent with 
observations (e.g. Hu et al. 1995; Kim ct al. 2002, and 
our own inferences, § 8.3). For the Doppler parameter 
distribution we set b a = 24 km s _1 (Kim ct al. 2001). 

Because SLLSs (19 < log(AT H i) < 20.3) and DLAs 
(log(A r Hi) > 20.3) have different column density dis- 
tributions and are usually excluded in measurements of 
T°g, these were added after the log(Ani) < 19 line for- 
est converged to the chosen t" s (z). To constrain the 
redshift evolution of SLLSs, we combined the sample 
by O'Meara et al. (2007) and the lower limit given in 
Rao et al. (2006), yielding £ S lls(z) ~ 0.066(1 + z) 1 - 70 . 
For the SLLS column density distribution we adopted 
/? = -1.4 (O'Meara et al. 2007). The DLAs were mod- 
eled via 4>la(-z) = 0.044(1 + z) 127 (Rao et al. 2006) and 
j3 = — 2 (Prochaska et al. 2005), ignoring deviations in 
f(Nui, X) from a single power law. The Doppler param- 
eter distribution was left unchanged. 

With the overall opacity of the modeled Lya forest con- 
sistent with observations, and the high column density 
systems incorporated, we used the generated line lists to 
compute H I Lyman series (up to Ly30) and Lyman con- 
tinuum absorption spectra. In total we computed 800 
different lines of sight for quasars in the redshift range of 
our sample, 160 each at z = 3.4, z — 3.6, z — 3.8, z = 4.0 
and z = 4.2, respectively. From these we then gener- 
ated mock SDSS spectra. First, the resolved H I forest 
spectra were multiplied onto synthetic quasar SEDs gen- 
erated from principal component spectra (Suzuki 2006). 
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Fig. 5. — This figure presents a sample of the mock spectra an- 
alyzed in the paper, (a): A spectrum representative of the typical 
data quality (S/N^ 12 ~ 2 at the Lyman limit of the quasar). There 
is a LLS (z = 3.88059, log N m = 18.27) that can be easily identified 
at this S/N^ 12 . Two systems at higher z with Nm ~ 10 16 cm~ 2 
modify the continuum but do not produce a T912 > 2 LLS; (b): A 
high-S/N spectrum (S/Ng^^ 6 at the Lyman limit of the quasar), 
rendering its partial LLS (2 = 3.98387, log A"ht = 17.24) easily vis- 
ible. An additional SLLS at lower redshift (2 = 3.76609, log Afjji = 
19.28) sets the flux to zero at the Lyman limit, (c): This mock 
spectrum shows a slow 'roll-off' in flux blucward of the quasar's 
Lyman limit. There is a LLS (2 = 3.68084, log A^ui = 17.87), and 
five systems with 16 < log A^hi < 17 at higher 2 that produce the 
roll-off. The S/N at the Lyman limit of the quasar is quite low 
(S/N~ 1.5) and the intervening systems further decrease the S/N, 
rendering the LLS invisible, (d): A spectrum with strongly under- 
subtracted sky background. The first strong system encountered 
is a LLS (2 = 3.87495, log Nm = 18.34), but the flux is above zero 
even after hitting a DLA (2 = 3.57678, log A^hi = 20.43). One can 
assess that the sky subtraction is poor, because the Lya profile of 
the DLA does not saturate. The solid/dotted veritcal lines in each 
panel trace the lyman limit/Lya line of each absorber. 

We then degraded the resolution of the mock spectra 
to R = 2000 by convolving them with a Gaussian, and 
rebinned them to Sv = 69km s , matching the approx- 
imate resolution and pixel size of the SDSS spectra. Fi- 
nally, we added Gaussian noise to the mock SDSS spec- 
tra. In each mock spectrum, the S/N was normalized in 
the quasar continuum at 1450A and varied as a function 
of flux and wavelength according to the throughput of the 
SDSS spectrograph. The sky level was approximated as 
a constant and readout noise was also incorporated. Fi- 
nally, for a subset of the mocks we imposed a sky subtrac- 
tion error implemented by subtracting/adding a constant 
to the spectrum. These were generated to mimic such 
systematic errors that occasionally occur in the SDSS 
spectra. 

Four representative examples are shown in Figure 5. 

7.2. LLS Recovery and Sky Subtraction Bias 

Two of the authors (JXP, JMO) analyzed mock spectra 
using the identical tools and procedures applied to the 
SDSS spectra; these steps were done without knowledge 
of the mock line distribution and column densities (con- 
structed by author GW). The integrated results for the 
two authors were nearly identical; the following discus- 
sion and figures refer to the results related to JXP. 

Figure 6 summarizes the completeness and several bi- 
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Fig. 6. — Offset in redshift Sz for the 'observed' LLS and 
pLLS from the z a b s value of the nearest 'true' T912 > 2 LLS in 
our mock spectra. The analysis is restricted to the highest 2 a b a 
LLS along each sightlinc. In the upper panel, we show the Sz 
value for sightlines where a LLS was 'observed' and actually exists. 
Lighter/darker (green/black) points correspond to LLSs discovered 
outside/ within the statistical survey path (i.e. S/N^ 12 =2). These 
are the majority of cases (> 80%) and we find small Sz values with 
a small, but important bias to Sz > 0.01. In the lower panel, we 
show false negatives (diamonds) and false positive (triangles) de- 
tections. The former correspond to a true LLS that was observed 
as (Sz < 0.1) or hidden by a pLLS (Sz > 0.1). The misidentifi- 
cations (Sz < 0.1) are dominated by the mock spectra with large 
underestimates in the sky background. The dominant effect of a 
sky subtraction bias is an underestimate in the incidence of LLSs. 
For the misidentifications with Sz < 0.1, the dark points indi- 
cate systems that satisfy all of the survey criteria (e.g. z > 3.3, 
S/NA a > 2). 

ascs uncovered by our analysis. In the top panel, all cases 
where a mock T912 > 2 LLS exists with z a b s > 3.2 and 
an LLS was 'observed' are presented; these correspond 
to > 80% ol the cases. Specifically, we plot the offset 
5z between the true and observed LLS absorption red- 
shifts as a function of the S/N^ 12 of the spectrum. We 
find excellent agreement (small Sz), nearly independent 
of the S/N^ 12 of the data. There are, however, a number 
of cases with 8z > 0.01, primarily related to the blending 
of absorption lines (see below) . 

The lower panel in Figure 6 presents those cases that 
are false negatives (diamonds) and false positives (trian- 
gles). The latter are very rare. The former, however, did 
occur quite often in our analysis and can be divided into 
two classes: (i) cases with \8z\ < 0.05 which are proper 
false negatives, i.e. true rgi 2 > 2 LLS that were modeled 
as a pLLS; (ii) cases with 5z > 0.1 which are sightlines 
where a higher z pLLS precluded the detection of a lower 
z LLS. The majority of the latter cases are due to bona- 
fide pLLS at higher z which greatly diminish the S/N 



of the spectra at shorter wavelengths and 'obscure' the 
presence of a Tgi 2 > 2 LLS. Almost none of these cases, 
however, satisfy the selection criteria established in § 6; 
either the data have too low S/N^ 12 or the absorption 
redshift of the LLS gives z cm — z a b s > 0.4. 

The first class of false negatives, meanwhile, are almost 
exclusively associated with spectra that had systemati- 
cally low estimates for the sky background. In these 
cases, a T912 > 2 LLS has an apparent flux at A < Agi2 
and therefore was modeled as a pLLS. This is the domi- 
nant effect of a sky subtraction bias. Although our mock 
spectra had even numbers of over and under-subtracted 
sky backgrounds, only the former are relatively easy to 
identify (large regions of spectra are significantly nega- 
tive) and ignore. The net effect of a random sky sub- 
traction error is a systematic underestimate in the inci- 
dence of LLS. We stress, however, that the magnitude 
and frequency of poor sky subtraction in the mock spec- 
tra was intentionally elevated so that we could explore 
these effects. The incidence of such effects within the 
SDSS spectra is much lower, an assertion supported by 
our follow-up spectra with Keck/LRIS (Appendix A). 
Therefore, we are confident that the sky subtraction bias 
gives rise to a less than 10% systematic error for t T >i{z). 

7.3. The pLLS Bias 

Originally, we intended to perform a search for LLSs 
in spectra with S/Ng 12 =l to maximize the pathlength of 
the survey (a nearly 4x increase over S/Ng\ 2 =2). Our 
tests with mock spectra and follow-up observations with 
Keck/LRIS (Appendix A) indicated that we could ro- 
bustly identify rgi 2 > 2 LLS in such data. We also noted, 
however, that many of the spectra showed pLLS candi- 
dates which reduced the S/N^ 12 to below 1 and made 
the search for rgi 2 > 2 LLS much more challenging. Our 
response was to redefine the search path by attenuat- 
ing the absorbed continuum due to any identified pLLS 
candidates and then reapply the S /N^ 12 criterion for the 
remaining z < z p lls spectral range. With this approach, 
the search path was frequently terminated by the pres- 
ence of a pLLS candidate. 

In principle, this modification should provide an unbi- 
ased search for T912 > 2 LLS. Our trials with mock spec- 
tra, however, revealed an insidious bias associated with 
this redefinition of the search path. Specifically, it is very 
difficult to identify pLLSs when they occur at a small 
redshift offset (Sz < 0.1) redward from a T912 > 2 LLS. 
In these cases, instead of the search being terminated at 
the redshift of the pLLS such that the lower z LLS is 
not included in the survey, it is the pLLS that is ignored 
and the LLS is included within the statistical sample. 
Furthermore, the redshift of the recovered LLS is biased 
to a higher value which reduces the survey pathlength 
by a small but non-negligible quantity. Together, these 
two effects lead to an overestimate of £ T >2(z) by values 
ranging from 30 — 50%. Furthermore, we find that it is 
very difficult to precisely estimate the magnitude of this 
systematic bias for it depends sensitively on l T >2{z), the 
incidence of pLLS, and the quasar z em distribution. In 
our opinion, one cannot robustly correct for this system- 
atic pLLS bias and we caution against performing any 
analysis that would be subject to it. For these reasons, 
we ignore pLLS when defining the survey path based on 
a signal to noise criterion and performing the search for 
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7912 _ 2 LLS. With this approach, one must have suf- 
ficient S/N^ 12 to identify LLS even when one or more 
pLLS modify the absorbed continuum. This last point 
motivated our decision to restrict the survey to spectra 
with S/N^ 12 > 2. 



7.4. The Blending Bias 



As noted in the previous sub-section, LLS and pLLS 
that lie close to one another in redshift are very difficult 
to distinguish as individual systems. This is even true 
in the limit where one has spectra with exquisite S/N 
and resolution when Sz < 0.1 (or less in the case of high 
resolution echelle observations). With SDSS spectra, the 
limited information provided by the Lyman limit and the 
strongest Lyman series lines is insufficient to robustly 
distinguish multiple LLSs from a single system. This 
leads to a "blending bias" that manifests itself in several 
ways. 

First, the blending bias increases the number of LLSs 
observed because pairs of pLLSs blend together to give 
a single system with T912 > 2. Second, the absorption 
redshifts of the LLSs are shifted to higher redshifts be- 
cause one generally adopts zlls from the higher of the 
pair of systems. This leads to a smaller survey path 
and possibly a higher inferred incidence of LLSs. More 
importantly (see below), many LLSs are shifted into the 
proximate region of the quasar. This causes an underesti- 
mate of £ T >2(z) for intervening LLSs and an overestimate 
of PLLSs." 

We explored the quantitative effects of the blending 
bias with the following analysis. We constructed a set 
of mock absorption lines for each quasar in the statisti- 
cal survey (Table 2) with an incidence set to match our 
measurements (§ 8.1). Specifically, we adopted an A^hi 
frequency distribution 



/(Afe, z) = CN^ 



1 



1 + z* 



0) 



119.5 , 



with f3 = -1.3 for 10 16 5 cm" 2 < 7V H i 
and j3 = -2 for N m > 10 19 5 , = 3.7, 7 = 5.1 and 
C = 1.244 x 10 5 . From the mock absorber list we iden- 
tified all LLSs with N m > 10 17 5 cm" 2 . This formed the 
control sample. Then, we blended together all systems 
with \Sz\ < 5zj where S Zj = [0.01,0.05,0.1,0.2] and rei- 
dentified systems satisfying Nm > 10 17 ' 5 cm~ 2 . When 
blending two or more systems together, we set z a hs to 
the maximum of all the lines. Finally, we calculated the 
incidence of LLSs using the survey path and LLSs for 
each Szj. 

Figure 7a presents the results of the blending bias in 
terms of the enhancement/decrement of the incidence of 
T912 _ 2 LLS relative to the perfect model (Szj = 0). 
For Szj < 0.1, there is only a small and ignoreable effect. 
For Szj > 0.1, however, we predict a systematic underes- 
timate for the incidence of intervening LLSs, especially 
at z > 4 where the absolute incidence is highest. This 
deficit in £ T >2( Z ) runs contrary to expectation and is en- 
tirely due to the redshift bias where the blended LLSs end 
up with zlls within the proximate region. In turn, we 
predict a systematic over-estimate of l r >i{z) for PLLSs 
(Figure 7b). 

Our experiments with mock spectra indicate that we 
commonly blend together systems with Szj — 0.05 to 0.1, 
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Fig. 7. — Blending bias (enhancement/decrement of i T >2( z ) 
relative to no bias) for mock absorption line systems for (upper) 
intervening T912 > 2 LLSs and (lower) proximate T912 > 2 LLSs. 
The absorption line statistics were set to roughly match the ob- 
served incidence. Even for blending with redshifts 8zj = 0.2, the 
bias for intervening systems is relatively small. In contrast, the 
blending bias systematically elevates £ T >2( Z ) for PLLSs, especially 
at z > 4. 

with a weak dependence on redshift or S/N. We conclude, 
therefore, that the effects of the blending bias on our 
SDSS analysis are < 10% for measurements of I T ~>i(z) for 
intervening LLSs. This is below the current level of the 
statistical error and we will ignore it in the presentation 
of the results and discussion. For the PLLSs, the effect 
for Szj =0.1 ranges from 20 — 90% and we will report 
our results on £ T >2(z) for these absorbers as upper limits, 
especially for z a b s > 4. 

Before closing this section, we stress that the blend- 
ing bias affects all previous and future LLS surveys. In 
particular, we caution that the incidence of LLSs cannot 
trivially be used to constrain the H I frequency distri- 
bution f(Nm,X), because the latter assumes that every 
absorption system is a unique, identifiable line. For the 
LLSs in particular (but this is also true for the Lya for- 
est) , blending smears these lines over non-negligible red- 
shift intervals (Sz « 0.1) and this affect must be consid- 
ered when comparing against theoretical line densities. 

7.5. Continuum Uncertainty 

An important systematic uncertainty in our analysis 
is the placement of the absorbed quasar continuum. As 
described in § 4, the continuum for each quasar spec- 
trum was determined from an automated fit of a template 
model to the data. Each continuum was then reviewed by 
two authors (JXP, JMO) and frequently scalled up/down 
by 5 — 10%. For our LLS survey, modifications to the con- 
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TABLE 4 
LLS STATISTICS 



AX Az a m£ LS z° i T > 2 {z) d Z T > 2 (z) d £ T > 2 (X) d Ar£ LS 

(S/NA 2 =2) (S/NA 2 =3) (S/NA 2 =2) (Mpc) 



Full Sample 



[3.30,4.40] 366.5 93.8 


192 


3.68 2.051°;^ 2.1lt° f s 


52+ ' 04 


78 


Redshift Intervals 



[3.30,3.40] 


25.7 


6.8 


9 


3.35 


, 01 +0.43 
1 - ol -0.60 


9 9 o+0.84 
z.zo_j 23 




U.OJ_Q 16 


147 


[3.40,3.50] 


49.2 


12.9 


28 


3.47 


9 17+0.40 
z - ± ' -0.49 


z.uo_ 73 


57+ - 11 
u -°'-0.13 


83 


[3.50,3.65] 


111.7 


29.0 


46 


3.58 


i.oy_ 27 


1 44+0-33 


41+0.06 
u - 41 _0 .07 


107 


[3.65,3.90] 


109.7 


27.9 


57 


3.74 


2 05+ ' 27 


9 O-+0.44 
*- ot -0.54 


52+ ' 07 
u - oz -0.08 


76 


[3.90,4.10] 


41.9 


10.4 


30 


3.97 


2.89±g:|| 


2.45±S;™ 


0.72+.°;" 


48 


[4.10,4.40] 


28.3 


6.8 


22 


4.23 


Q QO + 0.68 

J - zz -0.84 


3.44±i;§| 


0.78±°;|S 


38 



Note. — Unless specified otherwise, all quantities refer to the S/N^ 12 =2 threshold. The 
cosmology assumed has S7a = 0.7, H m = 0.3, and Ho = 72km s _1 Mpc -1 . 
a Total redshift survey path for the S/N^ 12 =2 criterion. 
b Number of T912 > 2 LLS discovered in the survey path. 
c Median absorption redshift of the LLS for the S/N^ 12 =2 threshold. 

d Line densities of LLS with r > 2 per redshift dz or absorption length dX. Often written as 
dn/dz (or dN/dX) in the literature. 
e Average proper distance between LLS with r > 2. 



tinuum primarily modify the survey path; the estimates 
of T912 are only affected if 7912 < 2. 

To test the sensitivity of our results to continuum 
placement, we reanalyzed the data after scaling each con- 
tinuum up/down by 5 and 10%. We also considered a sce- 
nario where the scaling was random between ±10%. We 
find that the estimates on the incidence of LLSs varies by 
5 — 10% as the continua are modified. The effect, while 
due to systematic changes in the continuum, is not sys- 
tematic. That is, a systematic increase in the continuum 
does not systematically increase/decrease the incidence 
of LLSs at all redshifts. Therefore, we conclude that con- 
tinuum placement errors yield a random, non-negligible 
(« 5 — 10%) uncertainty in the final results. 

7.6. Summary 

We have conducted an assessment of the systematic un- 
certainty related to surveying LLSs using mock spectra 
with idealized Lya forest absorption yet realistic spec- 
tral characteristics (noise, resolution). Our analysis re- 
vealed an insidious bias associated with pLLSs that is 
best minimized by restricting the analysis to data with 
S/Ng\ 2 > 2. We identified an unavoidable bias related to 
the blending of LLS and pLLS that implies a a « 10% 
uncertainty in the measured incidence of LLSs. This bias 
becomes even more significant at z > 4 when the inci- 
dence of LLSs exceeds 3 per unit Az. Finally, we found 
that continuum placement errors yield a random, non- 
negligible (« 5 — 10%) uncertainty. Although higher S/N 
and spectral resolution will reduce some of these effects, 
we conclude that it will be difficult to avoid a system- 
atic error of 10 — 20% using the standard approaches to 
surveying LLSs. We believe that future programs will 
require new techniques if higher precision measurements 
are desirable. 

8. RESULTS 



In this section, we present the principal results of our 
survey. We defer extended discussion of previous work 
and the implications of our analysis to the following sec- 
tion. Systematic biases and uncertainies in these results 
were discussed in the previous section and are summa- 
rized in § 7.6. 

8.1. £ T >2(z): The Incidence of Intervening 7912 > 2 
LLSs per Redshift Interval 

An LLS survey, by its nature, provides only a single 
observable quantity: the incidence of LLSs per redshift 
interval l T >i[z). This quantity is independent of any as- 
sumed cosmology and consequently has limited physical 
meaning. Nevertheless, it is the proper starting point for 
describing our results. 

Following standard practice, we estimate £ T >2(z) from 
the ratio of (-/Vlls) the number of LLSs detected in a 
redshift interval to (Aztot) the total search path for 
that redshift interval: 



Kt>2 



LLS 



(10) 



Aztot 

The statistical error in l T >i{z) from this estimator is as- 
sumed to be dominated by the Poisson uncertainty in 
•Wlls ■ We have discussed a range of possible systematic 
uncertainties with this estimator in the previous section. 
Figure 8 presents the values of £ T >2(z) for the S/Ng 12 =2 
criterion in a set of arbitrary redshift intervals chosen to 
give tolls ^ 30 systems per bin. Table 4 lists these val- 
ues for S/Ng 12 thresholds of 2 and 3; there is no obvious 
dependence with this threshold. 

Figure 8 reveals that the incidence of rg 12 > 2 LLS 
increases monotonically for z > 3.5. Following previous 
work, we have modeled the redshift evolution in £ T >2(z) 
as a power-law with the functional form: 



C-t>2(z) = Clls 



1 



1 + z* 



(11) 
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Fig. 8. — Incidence of intervening LLSs with rgi2 > 2 £ T >2( Z ) as a function of redshift (solid points). Only the darker points were 
included in a maximum likelihood analysis to determine best-fit power-law (blue curve): £ T >2(z) = Clls[(1 + z )/(l + z*)] 7LLS , with 
z* = 3.7, Clls = 1-9 i 0.2, and 7lls = 5.2 ± 1.5 (68% c.l.). The dotted points, meanwhile, show i T >2( z ) when one includes quasars with 
3.4 < z e m < 3.6. These measurements arc significantly biased to higher values by the SDSS quasar-targeting criteria (PWO09). In the 
subpanel, the black (solid) curve shows the cumulative number of T912 > 2 LLSs detected in our survey of the SDSS-DR7 database adopting 
the S/N^ 12 =2 threshold. The blue (dotted) curve shows the predicted number of LLSs assuming the best-fit power-law from Figure 8 
and adopting the g(z) function from Figure 4. A one-sided KS-test does not rule out the null hypothesis that the model distribution is 
statistically different from the observations. 

is our expectation that 7lls is likely lower than the cen- 
tral value of our analysis. This assertion will be tested 
with future observations at z < 3 and z > 4.5. 

For z < 3.5, Figure 8 shows two evaluations of £ T >2{z)- 
The light, solid points show the l r >i[z) values derived 
from our statistical quasar sample with the restriction 
that z em > 3.6. These values are consistent with an ex- 
trapolation of the best-fit power-law. The dotted points 
in the figure, meanwhile, show the values of £ T >2{z) when 
one also surveys quasars with 3.4 < z cm < 3.6. In this 
case, we find systematically higher £ T >2(z) values which 
would indicate a non-physical, non-monotonic evolution 
in l T >i{z). These results confirm the findings of PWO09 
that the SDSS targeting criteria for quasar spectroscopy 
biases the sample against sightlines without a LLS. The 
values of £ T >2{z) reported in Table 4, therefore, are re- 
stricted to quasars with z em > 3.6. 

8.2. The Incidence of LLSs in ACDM 

If one introduces a cosmological model, the observed 
incidence of 7912 > 2 LLSs with redshift £ T >2{z) may 
be translated into physical quantities. Consider first, 
ArLLS) the average distance that a photon travels be- 



setting z* = 3.7. Using standard maximum likelihood 
techniques (e.g. Storrie-Lombardi et al. 1994), we find 
best-fit values to the data at z > 3.5 of Clls = 1-9 ± 0.2 
and 7 L ls = 5.2 ± 1.5 (68% c.l.). The best-fit model is 
overplotted on the data in Figure 8. The relatively large 
uncertainty in 7lls is due to the small redshift interval 
covered by our survey. Nevertheless, we conclude at high 
confidence (> 95%) that £ T >2(z) is increasing at least as 
steeply as 7lls = 2 at z > 3.5. 

The sub-panel of Figure 8 compares the cumulative 
number of LLSs detected in the survey against red- 
shift both as observed (solid) and as predicted (dot- 
ted) by the best-fit power-law model. For the latter, 
we adopt the g(z) curves for S/N^ 12 =2 from Figure 4. A 
one-sided Kolmogorov-Smirnov test yields a probability 
Pks = 0.95 that the observed distribution is drawn from 
the adopted power-law expression; the power-law model 
is a good description of the observations. We comment, 
however, that the best-fit slope (7lls = 5-2) is consider- 
ably steeper than most previous estimates for the LLSs 
at this redshift (see § 9.1) and also steeper than the red- 
shift evolution measured for the Lya forest and damped 
Lya systems (Kim et al. 2001; Prochaska et al. 2008). It 
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Fig. 9. — Incidence of rgi2 > 2 LLSs as a function of absorption path X as a function of redshift, indicated by the dashed lines. We 
measure a decrease of 2 in ^ T >2p0 from z = 4.4 to 3.5. The widths of the two lower bands (blue, red) indicate estimates of £(X) for 
DLAs and SLLSs (O'Meara et al. 2007; Prochaska & Wolfe 2009); see the text for details. These bands are plotted on top of one another 
to indicate their total contribution to t T ->i{X) (as indicated by the vertical arrows on the right-hand side of the figure). The (top) blue 
band, therefore, represents the estimate to l T >2(X) from LLSs with Nm = 10 17 ' 5 — 10 19 cm - 2 . Given the low absolute value and weak 
evolution in £(X) for the SLLSs and DLAs, the evolution in £(X) for the LLSs is likely dominated by systems with T912 < 10. 



fore encountering an LLS with T912 > 2. Specifically, we 
define 



where 



and 



Ar LLS =£ T>2 (z)- 1 ^ , 
az 



dr 



dz (l + z)H(z) 



h{z) = h q [n A + (1 + zfn m ] 



1/2 



(12) 
(13) 
(14) 



With our adopted cosmology, we estimate that Arrxs 
ranges from « 100 to 40 Mpc proper distance from 
z = 3.5 to z = 4.4 (Table 4). This is an order of 
magnitude or more larger than the separation of high 
z quasars (several Mpc for Lb > 10 40 ergs -1 Hz -1 ; 
Faucher-Giguere et al. 2009). 

An especially informative quantity for associating LLSs 
to structures in the Universe (e.g. galaxies, filaments) is 
£{X) the number of systems per absorption length (Bah- 
call & Peebles 1969), where £(X)dX = £{z)dz and 



?dz 



(15) 



The quantity £ T >2(X) is defined to remain constant if 
"■lls 7 the comoving number density of structures giving 
rise to a rg 12 > 2 LLS, times ullSi the average physi- 
cal size of the structure remains constant, i.e. £ T >2(X) oc 
"llsclls- Figure 9 presents the evolution of £ T >2(X) for 
our cosmology as a function of redshift (see also Table 4) . 
We observe a rise in £ T >2(X) with redshift of roughly two 
times over the rj 1 Gyr from z = 3.3 to 4.4. At 99% con- 
fidence, we infer an increase in £ T >2{X) over this redshift 
interval. This follows, of course, from the very steep red- 
shift evolution observed for £ T >2{z) (§ 8.1); in a flat cos- 
mology with fl nl on the order of Oa, an £ T >2{z) evolution 
steeper than (1 + z) 1 / 2 implies £ T >2(X) is also increasing. 
We conclude that tt-lls and/or ctlls are increasing with 
redshift at z w 3.5. We discuss the implications of this 
result in S 9.4. 



8.3. f(N m ,X) at z « 3.7 

In this subsection, we combine our results with pre- 
vious work on the IGM to place constraints on the H I 
frequency distribution, f(Nm,X). We focus this analy- 
sis at a single redshift (z = 3.7) where our observations 
have greatest statistical power. 
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Fig. 10. — The Nm frequency distribution f(Nm,X) observed for the SLLSs (green, iVjji 



10 1 



" 2 : O'Meara et al. 



2007) and DLAs (red, 
N m = 10 17 ' 5 - 10 19 cm 



^Vhi > W cm' 2 ; Prochaska & Wolfe 2009). The blue band is an estimate of f(N m ,X) for LLSs having 
~ 2 under the assumptions of a power-law form (f(Nm,X) oc Ar£i' :LS ) and constrained by the observed incidence of 



SLLSs and r > 2 LLS (this paper). We find /3 L ls 



-0.8 



+0.4 



(68% c.l.) for conservative estimates on the value of f(Nm, X) at 10 cm 



and allowing for 20% uncertainty in £ T >2{X). The dashed and dotted curves indicate the range of power-laws that satisfy the observations. 



8.3.1. f hhS (N m ,X) 

Although our observations and LLS analysis are in- 
sensitive to the H I column densities of the LLSs, they 
do provide an integral constraint on the frequency dis- 
tribution of Nbi, f(Nni,X). We constrain the Am 
frequency distribution of T912 < 10 LLS per absorp- 
tion length, /lls(Ahi,A), at column densities Ahi = 
lO 17 5 — 10 19 cm~ 2 as follows. Previous surveys at 
z > 3 have measured /(Nm, A) for column densities 
A H i > 10 19 citT 2 (Prochaska et al. 2005; O'Meara et al. 
2007; Prochaska & Wolfe 2009; Noterdaeme et al. 2009; 
Guimaraes et al. 2009). These authors have parame- 
terized the distribution functions as single (SLLSs) and 
double (DLAs) power-laws of the following form: 



/ S LLs(10 19 cm- 2 < N m < 10 20 - 3 cm- 2 ,A) 



k 



;slls<i lls 



and 



(16) 



/dla(A H i > 10 2, ^cm- 2 ,A) = fc DLA 



Nm 



where /3dla 



/3 3 : A H i < N d ; 
04. ■ Ami > N d 



(17) 



Figure 10 presents these frequency distributions. For the 
SLLSs at z = 3.7, we have taken /?slls = —1.2 ±0.2 and 
normalized the power-law by taking 



^SLLS(A) 



10 1 



/sLLs(A H i,X)dA/ H i = 0.2 



(18) 

These values are consistent with the range of pub- 
lished measurements at this redshift (Peroux et al. 2005; 
O'Meara et al. 2007; Guimaraes et al. 2009). For the 
DLAs, we have evaluated /dla(Ahi, X) from the SDSS- 
DR5 (Prochaska & Wolfe 2009) over the redshift interval 
z = [3.4,4.0], giving N d = 10 21 - 75 , ft = -1.8, ft < -3, 
and /cdla = 7 x 10 25 cm 2 . 



f(Nm,X) for the interval A ffl = 
' which we refer to as /lls(Ahi, A), 



We estimate 
[10 17 - 5 , 10 19 ] cm 
under the following assumptions/constraints 
/lls(Ahi, A) has a power-law form 



/ L Ls(A r Hi,A) = fc LLS A ] 



HI 

10 19 ' 



(19) 



and (ii) /lls(Ahi,A) at Nya = 10 ia cm z is consistent 
with the range of values given by the SLLSs. Specifically, 
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we demand log/(A ffl = 10 19 cm- 2 ,A) = -20.05 ± 0.2; 
(iii) we impose the integral constraint based on the the 
observed incidence of Tgi 2 > 2 LLSs: 



ir>2(X) = 



f(N m ,X)dN ] 



HI 



(20) 



in 1 



At z k. 3.7, we estimate l r >i(X) = 0.5 ± 0.1 (Figure 9). 

Overplotted on Figure 10 are the power-law frequency 
distributions (shown as dashed, dotted and a solid line) 
that satisfy the extrema of those constraints. The shaded 
region shows the intersection of the curves and roughly 
represents the allowed region of /(Am, A) values. We 



10 



-4.5 



cm 



find /3lls = —0.8 ± 0.3, and derive /clls 
for the central value. Table 5 further summarizes these 
results. 

Our analysis reveals that /(Am, A) becomes increas- 
ingly shallow with decreasing Am- Only for the most 
extreme values of our analysis, low /(Am, A) at Am = 
10 19 cm -2 and a large £ T >2{z) value, do we recover 
/?lls < — 1- This flattening of /(-/Vhi, A") was sug- 
gested by previous authors based on a similar analy- 
sis but with much poorer observational constraints on 
£ T > 2 (X) (Peroux et al. 2003; Prochaska et al. 2005; 
O'Meara et al. 2007). Remarkably, our results indicate 
/?lls > — 1 which means that the IGM has a higher total 
covering fraction per unit pathlength for sightlines with 
N m = 10 19 cm -2 than those with A H i 



10 18 cm" 2 . 



8.3.2. Constraints from Measurements of the Mean Free 

Path 

Traditionally, the mean free path to ionizing radiation 
in the IGM (A 9 ^) has been estimated from the observed 
incidence of LLSs (e.g. Meiksin & Madau 1993). Re- 
cently, PWO09 introduced a new approach to measure 
^mfp fr° m stacked quasar spectra, without any considera- 
tion of LLSs. We can reverse the problem, therefore, and 
use the A 9 ^ results to constrain properties of /(Am, A). 
One expects to have the greatest sensitivity to absorption 
systems with r 9 i 2 ~ 1, i.e. the LLSs and partial LLSs. 

At z = 3.7, PWO09 estimate A^ 2 = 47/if 2 1 Mpc 
proper distance. This means that in the absence of an ex- 
panding universe, a packet of 1 Ryd photons at z = 3.7 
would be attenuated by exp(— 1) after traveling X^ p . 
One can also express the mean free path as an opacity, 
K 9i2=l/A 9 n 1 f 2 , which can be related to the optical depth 
of a 1 Ryd photon as: 



0.025 



0.020 



X 

Z 
A 



0.015 



a, 
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Fig. 11. — This figure shows the opacity at the Lyman limit 
.12 contributed by absorbers with N m > N™ n . For N™ n > 
17.5 cm -2 ( so ii(j curve), which corresponds to our LLS survey, 



we have adopted the estimate of f(Nn-[,X) from Figure 10 in 
the calculation. We estimate a 20% uncertainty in the contri- 
bution of LLSs to K912, as shown in the figure. For iV5j in < 
10 17,5 cm -2 , we assume f(Nm,X) follows a simple power-law 
with exponent /3 p llS; an d show a series of extrapolations (dash 
and dotted curves). The solid (red) horizontal band centered at 
K912 = 0.0225 hf2 Mpc -1 indicates the measurement at 2 Ri 3.7 by 
PWO09. These results imply that LLSs contribute rj 55% of the 
mean free path (^ 33% for systems with rgi2 ^> 1) and that /3 p lls 
must be steeper than « —1.5 to explain all of these observations. 

G?r 9 jP is insensitive to the minimum Ahi column density 
for any value Af : in < 10 12 cm" 2 . 

Figure 11 shows the K912 value at z = 3.7 from 
PWO09 as a horizontal band that illustrates the la er- 
ror interval. The solid curve, meanwhile, corresponds 
to the evaluation of Equation 21 using our best estima- 
tion of /(Ahi, A) (Figure 10) as a cumulative function 
of A^j 111 . At the limiting Ahi value of our LLS survey 
(10 17 5 cm -2 ), we estimate that w 55% of the opacity 
to ionizing radiation is contributed by r 912 > 2 LLSs. 
The uncertainty in the results is roughly proportional to 
the uncertainty in £ T > 2 {z), i.e. ss 20% as indicated by 
the error bars on the figure. It is notable that « 1/3 
of the contribution to Kgi 2 is from very optically thick 
absorbers (rgi 2 ^ 1), i.e. the SLLSs and DLAs. 

It is also evident from Figure 11 that systems with 
T 9i2 < 2 must contribute to Kgi 2 . For A2j n < 



10 cm , we continue the calculation by assuming 
that /(Ahi, A) follows a power-law 



K912 



dr^ _dr^dz 



dr 



dz dr 



(21) 



Finally, we can relate the differential optical depth to the 
H I frequency distribution of absorbers: 



dz 



00 

J /(Am, z){\ - cxp [ -Ahict 9 ^ ] }dA H i 



(22) 

with cr 9 ^ 2 the photoionization cross-section evaluated at 
lRyd and dr j dz given by Equation 13. Although the 
integral should be evaluated with Agj ln — 0, in practice 



/pLLS (Ahi < 10^ ^cm- 2 ,A) = fe pLLS A^ 



ftjLLS 



(23) 



constrained to match /lls(Ahi,A) at Ahi — 
10 17 - 5 cm- 2 . Wc find that models with /3 pLLS > -1.8 
cannot reproduce the A 9 ^ 2 results. In fact, the data favor 
/3pLLS ~ — 2 i.e. a much steeper power-law than inferred 
for the LLS and also than that commonly observed for 
the Lya forest. These conclusions depend rather insensi- 
tively on our estimate of £ T > 2 (A); slopes only as shallow 
as —1.7 are allowed if we adopt our highest estimates for 

ir>2(X). 

Thus far, these inferences on /(Ahi, A) for absorp- 
tion systems with r 9 i 2 < 2 have ignored observations 
of the Lya forest. By including these data, we provide 
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Fig. 12. — The solid black curve shows our estimation of f(Nm,X) at z Kb 3.7 as a series of six power-laws that intersect at Afjji = 
[10 14 5 , AfpLLs, 10 19 , 10 20 3 , 10 2175 ] cm" 2 , where 7V pLLS is constrained to lie between jV HI = 10 15 - 10 17 5 cm- 2 . The observational 
constraints for f(Nm,X) at Nm > 10 19 cm -2 are as in Figure 10. The results here also include constraints from the observed mean free 
path (PWO09) and the effective Lya opacity of the Lya forest (Faucher-Giguere ct al. 2008b). The orange curves show the regions of 
f(Nm,X) that reproduce the mean free path measurement and also connect to the Lya forest at TVhi = 10 14,5 cm -2 . 



further constraints on /(Am, X) for Am — 10 15 cm~ 2 
to 10 19 cm -2 . To derive these constraints, however, we 
must adopt a functional form for /(Am, A). Absent a 
physical model, we take an empirical approach. We ex- 
press /(Am, A) as a series of six power-laws that inter- 
sect at iV ffl = [10 14 - 5 , iV pLL s, 10 19 , 10 20 - 3 , 10 21 - 75 ] cm" 2 , 
where -/V p lls is constrained to lie between Am = 10 15 — 
10 17 ' 5 cm -2 . Other than a small 'kink' at Am = 
10 20 ' 3 cm" 2 , the power-laws are required to match at 
each intersection point. 

The power laws are forced to satisfy the following ob- 
servational constraints: 

1. The power-laws for Am > 10 175 are constrained 
as described at the start of this sub-section. 

2. For the Lya forest, we assume /(Am, A) oc A ffl 15 
and normalize at z = 3.7 by the effective Lya opti- 
cal depth t" s measured by Faucher-Giguere et al. 
(2008b) assuming a &-value distribution f(b) oc 
&- 5 exp[-& 4 /6 4 ] ( Hui & R u tledge 1999). 

3. The integrated opacity of the IGM at the Ly- 
man limit is constrained by the measurement of 
PWO09, i.e. = 47 Mpc. 



For the range of power-laws derived from our analysis 
of the LLS results (Figure 10), we show in Figure 12 
the range of /(Am, X) distributions that also satisfy all 
of the constraints. We find viable models with A p lls 
values that range across the allowed interval. These are 
correlated with /3 p lls values ranging from /3 p lls = — 1-9 
to -5 (Table 5). 

The principal results of this analysis are threefold. 
First, the single power-law connecting the Lya for- 
est to the SLLS satisfies neither the LLS nor mean 
free path constraints. There is at least one break be- 
tween Am = 10 14 - 5 cm -2 and Am = 10 19 cm~ 2 where 
/(Am, A) steepens to f3 < —1.8 and then flattens to 
/9w-l. This is consistent with conclusions drawn from 
line-counting statistics of Lya forest lines (e.g. Petitjean 
et al. 1993; Kim et al. 2002). Second, the added A 9 ^ 2 
and Lya forest constraints rule out the lowest values of 
/lls(Ahi,A) at A H i ~ 10 17 - 5 cm -2 that were otherwise 
allowed by our LLS results. Specficially, the data re- 
quire log /(A H i = 10 17 cm~ 2 ,A) > -19 and the (blue) 
shaded region in Figure 12 shows the proper allowed 
range for /lls(Ahi,A). Finally, we find that the slope 
of /(Am, A) must steepen to /3 p lls < —1-8 at columns 
Am > 10 14 - 5 cm~ 2 . 
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TABLE 5 
SUMMARY OF f(N m ,X) RESULTS 



«LLSv*)° log/f 9 ft,LS logfcLLS log^pLLs" /3pLLS 



Preferred Values 



0.23 -20.05 -0.8 -4.5 17.3 -1.9 



TABLE 6 
SDSS-DR7 PROXIMATE rgi 2 > 2 
LLS SURVEY 

Quasar z om z LL s 

J001115. 23+144601. 8 4.967 

J001820.71+141851.5 3.936 

J004219.74-102009.4 3.880 

J004240.65+141529.6 3.687 3.684 

J010619.24+004823.3 4.449 

J011351.96-093551.0 3.668 

J012403. 77+004432.7 3.834 

J014049. 18-083942.5 3.713 3.693 

J015048.82+004126.2 3.702 

J015339.61-001104.8 4.194 

J021318.98-090458.3 3.794 3.797 

J022518.35-001332.2 3.628 

J024447.78-081606.1 4.068 

J025518.58+004847.6 3.989 

J031213.97-062658.8 4.031 

J034402.85-065300.6 3.957 

J073149.50+285448.6 3.676 

J074154.59+341252.1 3.905 

J074500.47+341731.1 3.713 

J074640.16+344624.7 4.010 

J074711. 15+273903.3 4.154 

J075006.62+491834.1 3.603 

J075017.49+405825.3 3.864 3.849 

J075103.95+424211.6 4.163 

J075347.41+281805.2 4.031 

J075552.41+134551.1 3.673 

J075732.89+441424.6 4.170 



Conservative Range of Allowed Values 



0.23 


—20.05 


—0.8 


—4.5 


17.1 


—2.0 










17.3 


-1.9 










17.5 


-1.9 


0.15 


-20.25 


-0.9 


-4.1 


17.5 


-1.9 


0.35 


-20.25 


-1.3 


4.7 


15.0 


-3.5 










15.2 


-3.0 










15.4 


-2.6 


0.15 


-19.85 


-0.1 


-18.0 






0.35 


-19.85 


-0.8 


-4.6 


15.2 


-5.2 










15.4 


-4.3 










15.5 


-3.7 










15.7 


-3.3 










15.9 


-3.0 










16.1 


-2.7 










16.2 


-2.5 










16.4 


-2.4 










16.6 


-2.2 










16.8 


-2.1 










17.0 


-2.0 










17.1 


-1.9 










17.3 


-1.9 



Note. — The analysis throughout assumes that the DLAs 
contribute t(X) = 0.09 to £ T > 2 (X). 



a The incidence of T912 > 2 LLS t T >2(X), attributed to LLS 
with JVhi < 10 19 cm -2 . 

b The adopted value of f(N m ,X) at N m = 10 19 cm- 2 . 

c The break column density within the Lya forest as defined in 
the text. Entries without values have power-law descriptions 
for the LLS that cannot satisfy the mean free path and Lya 
forest constraints. 

8.4. The Incidence of Proximate LLSs (PLLSs) 

The results presented thus far all refer to intervening 
LLSs, i.e., systems restricted to have zlls blueward of 
3000km s _1 from the quasar emission redshift. This re- 
striction was imposed to isolate the 'ambient' IGM and 
avoid biases related to having performed the search for 
bright, background quasars. In the space surrounding 
a bright quasar, one predicts at least two such biases: 
(1) bright, high z quasars are known to cluster strongly 
(r > 15/i _1 Mpc; Shen et al. 2007) suggesting these 
objects trace massive structures in the young universe. 
The local environment of bright quasars, therefore, has 
an uncommonly high density (at least in dark matter) 
which may give a higher incidence of LLSs; (2) the ra- 
diation field of the quasar will ionize gas to large dis- 
tances, reducing the incidence of LLSs. For damped Lya 
systems, the first effect dominates as one observes an en- 
hanced rate of proximate DLAs (PDLAs) relative to the 
intervening systems (Russell et al. 2006; Prochaska et al. 
2008). 

Following the formalism presented in Prochaska et al. 
(2008), we have estimated the incidence of PLLSs in a 
series of redshift intervals. First, we re-measured the 
quasar emission redshifts for all systems with S/N^ 12 > 2 
at the Lyman limit and with a T912 > 2 LLS within 
5000km s _1 of z om . The SDSS quasar redshifts reported 
in the standard DR7 data release are known to have sig- 



Note. — [The complete version of 
this table is in the electronic edition of 
the Journal. The printed edition con- 
tains only a sample] 

nificant systematic errors. Following the prescriptions 
described in Shen et al. (2007), J. Hcnnawi has kindly 
remeasured the redshifts for all of the quasars. Sec- 
ond, we analyzed the quasars whose absorbed contin- 
uum at the wavelength of the Lyman limit correspond- 
ing to 3000km s _1 blueward of z om is twice the median- 
smoothed, la error-array. This establishes the survey 
path. All PLLSs identified redward of this 3000km s _1 
offset form the statistical sample (Table 6). The inci- 
dence, ^plls(z), is then estimated in arbitrary redshift 
intervals assuming the same estimator for intervening 
LLSs (Equation 10). These results are presented in Fig- 
ure 13 and compared against the incidence of intervening 
LLSs. 

Ignoring the data at z < 3.6 (which we suspect to 
be biased high by the SDSS targeting criteria; PWO09), 
the incidence of PLLSs roughly tracks that of intervening 
LLSs but is ps 25% lower than the intervening systems. 
The inset figure shows the observed cumulative number 
of PLLSs (dark curve) versus the predicted number (light 
curve) assuming the best- fit power-law for £ T >2{z) of in- 
tervening LLSs. A one-sided KS test yields only a 1% 
probability that the two distributions are drawn from 
the same parent population. We also remind the reader 
that corrections for the blending bias described in § 7.4 
will likely reduce £plls( z ) further, especially at z <~ 4. 
The principal implication is that T912 > 2 LLSs toward 
i < 20 mag quasars at z > 3.5 suffer from a proxim- 
ity effect, presumably due to the ionizing radiation field 
of the quasar itself. Given the observed enhancement 
of strong LLSs at z ss z em along sightlines transverse 
to such quasars, our results lend further evidence that 
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Fig. 13. — The solid and dark points show the incidence ^pllsC- 2 ) 
of proximate LLSs (PLLSs; LLSs with rgi2 > 2 that occur within 
3000 km s~ x of the background quasar) per unit redshift against 
quasar emission redshift. These are compared against the same 
quantity for intervening LLSs (gray points). The data point at 
2cm < 3.6 has been dotted out because it is biased by the SDSS tar- 
geting criteria (PWO09). Ignoring that last point, we find that the 
incidence of PLLSs roughly follows that of intervening systems but 
is systematically lower by fa 25% at z < 4. The inset figure shows 
the cumulative number of PLLSs observed (black curve) against the 
predicted number using the power-law model for £ T >2{ Z ) (§ 8-1) 
and the g(z) curve for PLLSs (Figure 4). 

quasar emission at « 1 Ryd is usually anisotropic (Hen- 
nawi & Prochaska 2007). 

Before concluding this section, we comment that the 
PLLS analysis is subject to another systematic error. In 
performing our LLS survey, we have identified and re- 
moved all quasars with very strong associated absorp- 
tion, e.g. BALs. The intent of this procedure was to 
remove the signatures of absorption from gas very local 
(< lkpc) to the quasar from the analysis. It is possible, 
however, that the associated absorption in some of these 
removed quasars is due to a PLLS at distances 3> 1 kpc 
and not very local gas. This would lead to an under- 
estimate of ^plls(z)- Alternatively, we may not have 
identified all of the local absorbers and therefore might 
have overestimated ^plls^)- I n either case, we caution 
that a systematic error of the order of 10 — 20% should 
be attributed to this effect. 

9. DISCUSSION 

9.1. Comparisons with Previous Work 

Surveys for Lyman limit systems have been carried out 
for several decades now (Tytler 1982; Sargent et al. 1989; 
Lanzetta 1991; Storrie-Lombardi et al. 1994; Stengler- 
Larrea et al. 1995; Peroux et al. 2003). These have been 
performed primarily at optical wavelengths on heteroge- 
neous quasar samples drawn from a diverse set of sur- 
vey approaches: color-selection, radio detection, slitless 
spectroscopy, etc. The authors adopted differing com- 
pleteness limits for T912 (ranging from 1 to 3) and used 
different approaches to establishing the pathlength that 
establishes l T >i[z). Little attention was given to assess- 
ing systematic error, and several of the effects described 
in § 7 assuredly apply to the previous works. Perhaps not 
surprisingly, therefore, most of the previous estimates of 
l(z) are in disagreement with our results. 



Fig. 14. — Estimates of £ T >i(z) from several previous stud- 
ies, parameterized as power-laws with the parameters listed in Ta- 
ble 7. The dotted (gray) lines show estimates from Sargent et al. 
(1989) and Lanzetta (1991). The dashed curve is the estimate from 
Stengler-Larrea et al. (1995) and the solid blue and green curves 
are from Storrie-Lombardi et al. (1994) and Peroux et al. (2003) 
respectively. For the values of our survey, we have increased the 
^t>2( z ) results by 10% to match the T912 = 1 threshold of the 
previous work. Our results indicate a significantly lower incidence 
of LLSs at z < 4 than suggested by the previous estimates. This 
lower incidence is consistent with recent estimates of the mean free 
path to ionizing radiation (POW09). 



TABLE 7 

POWER-LAW PARAMETERS OF i T >i(z) 



7LLS Clls Reference 

0.68 0.76 Sargent et al. (1989) 

5.7 0.00081 Lanzetta (1991) 

2.8 0.043 Storrie-Lombardi et al. (1994) 
1.5 0.25 Stengler-Larrea et al. (1995) 
2.45 0.07 Peroux et al. (2003) 



In Figure 14 we present estimates of l(z) from sev- 
eral previous studies, parameterized as in Equation 11 
(Table 7 lists the parameters) . The dotted lines show es- 
timates from Sargent et al. (1989) and Lanzetta (1991), 
which included very few observations at z > 3.5 and 
had very discrepant guesses for the high z universe. The 
dashed line shows the result from Stengler-Larrea et al. 
(1995) who integrated previous work with a new mea- 
surement at z < 1 and an (still) unpublished survey by 
Steidel & Sargent. Finally, the solid curves show the 
results from Storrie-Lombardi et al. (1994) and Peroux 
et al. (2003) who surveyed LLSs at z ~ 4 using color- 
selected quasars. All of these analyses were claimed 
to correspond to the incidence of LLSs with T912 > 1, 
^t>i(z), although a careful review of the literature raises 
doubts regarding this assertion. Nevertheless, to make 
comparisons with their reported T912 > 1 results we have 
boosted each of our l T >i{z) estimations. Formally, we 
estimate a correction of 6% from our derived f(Nm,X) 
distribution but, in practice, we adopt a more conserva- 
tive 10% correction. 

Our results on £ T >i(z) at z > 4 are lower than the 
values derived from the APM surveys, but within w la 
of concordance. The more important differences are be- 
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tween the estimations at z < 4. All of the previous work 
was essentially derived from the surveys of Sargent et al. 
(1989) and Lanzetta (1991) and, therefore, the curves all 
intersect at z « 3 at a t r >\(z) ~ 2. Our results suggest 
that much of the previous work at z ~ 3 overestimated 
the incidence of LLSs. The original survey by Sargent 
et al. (1989) is in fair agreement with an extrapolation 
of our power-law form for i T >x{z) to z = 3, but the re- 
analysis by Stengler-Larrea et al. (1995) of unpublished 
spectra taken by Steidel & Sargent led to a higher esti- 
mate at this redshift. We suspect that this later work 
gave t T >\{z) values that were too high either because of 
sample variance (i.e. small number statistics), selection 
bias in the quasar sample, and/or the systematic effects 
described in § 7. 

It is also reasonable to consider whether our SDSS sur- 
vey has been biased low by an unidentified systematic 
error. The measurement of A^ f p, however, indicates that 
this is not the case. Reconsider the analysis presented 
in § 8.3.2 (Figure 11). If we adopted a 50% higher inci- 
dence of T912 > 2 LLSs, e.g. £ T >2(X) — 0.75 at z = 3.7, 
then we would infer a much steeper slope for the LLSs 
(Plls ~ —1-5) and then would require a power-law shal- 
lower than /3 = -0.5 for N m = 10 15 cm- 2 - 10 17 cm- 2 . 
This would force f(Nm,X) to steepen to (3 < — 3 at 
TVhi ~ 10 14 - 5 cuT 2 '. Such an extreme f(N m ,X) distri- 
bution is non-physical and, more importantly, ruled out 
by observed line-statistics of the Lya forest (Kim et al. 
2002; Misawa et al. 2007). We conclude that the in- 
cidence of LLSs at z » 3.7 cannot be more than 30% 
higher than our central value, and, at present, we cannot 
identify a bias that would lead to such a large systematic 
underestimate in our results. 

9.2. Evolution in £ T >i(z) from 2 = — 4 

Previous work has debated whether l T >\{z) evolves as 
a singe power-law (l+z) 7LLS from z » 0—4 (e.g. Stengler- 
Larrea et al. 1995). In § 8.1, we modeled our observations 
with a single power-law having 7lls = 5.2 ± 1.5. The 
observed evolution in the mean free path (PWO09) also 
suggests a steep evolution (7lls > 2) for the LLSs. We 
now consider whether a single power-law extrapolation 
is a good description of £ T >x(z) for z < 3. 

Survey of other H I absorption systems have demon- 
strated that a single (1 + z) 7 power-law is a poor de- 
scription of the Lya forest (Weymann et al. 1998) and 
the damped Lya systems (Prochaska et al. 2008). In the 
former case, one observes a flattening in the Lya line- 
density at z ~ 1 which has been interpreted to result 
from a corresponding decline in the intensity of the extra- 
galactic UV background (EUVB Weymann et al. 1998; 
Dave et al. 1999). It is plausible that a similar effect 
would influence the LLS. In Figure 15, we present esti- 
mates of £ T >i(z) at z ~ 3.5, estimated by increasing the 
measured £ T >2(z) values by 10%. These are compared 
against the z ~ 1 measurement from (Stengler-Larrea 
et al. 1995, see Ribauldo et al., in prep for a new es- 
timate). Overplotted on the data is a solid curve that 
shows the best-fit power-law to the SDSS results. The 
dashed curves show 2a departures from this model, ex- 
trapolated to z = 0. It is evident that none of these 
curves intersect the low redshift observations. We con- 
clude, with high confidence, that a strict (1 + z) 7LLS 
power-law does not describe the evolution of £ T >i(z) from 
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Fig. 15. — A comparison of our results from the SDSS-DR7 
survey of LLSs (scaled to T912 > 1) against the z < 1 results from 
Stengler-Larrea et al. (1995). The (blue) solid and dashed curves 
show the best-fit to the SDSS observations and ±2c deviations 
from that fit. None of the curves, when extrapolated to z = 0, 
intersect the low z observations. We rule out at high confidence 
that a single power-law (1 + z)tlls ; s a good description of the 
observations from z = — 4. Instead, we suggest a break at z ^ 2, 
here modeled (dotted, red line) as two power-laws with 7lls = 2.78 
for z > 2.3 and 7lls = 0.26 otherwise. 

z = — 4. Instead, the data suggest a 'break' in the high 
z power-law at z ~ 2, similar to that observed for the 
Lya forest although at somewhat higher redshift. 

The dotted curve, is an attempt to model this break. 
For z < 2.3, we adopt the power-law form that matches 
the Lya forest at low redshift (7lls = 0.26) and demand 
that it intersect the central value of the Stengler-Larrea 
et al. (1995) measurement. For z > 2.3, the model breaks 
to a 7lls — 2.78 power-law, again consistent with the 
high z evolution of the Lya forest (e.g. Kim et al. 2002). 
This is a reasonably good description of the data and 
we conclude that a break in the power-law description 
of i T >i{z) likely occurs at z sa 2. We will test this pre- 
diction with an (ongoing) survey for LLSs at z ~ 2 in 
HST/ACS and WFC3 slitless spectra of z « 2.3 quasars 
(PI: O'Meara). 

9.3. The Average LLS Spectrum 

To gain additional insight into the absorption prop- 
erties of the LLSs, we have constructed an average 
(stacked) spectrum by (i) shifting each quasar spectrum 
containing a T912 > 2 LLSs to its rest-frame (192 systems 
total) and (ii) averaging the fluxed data. This stacked 
spectrum is primarily illustrative; it is shown in Fig- 
ure 16. The peak in emission at A « 1280A is from the 
Lya emission peak of the background quasars. The peak 
is offset from 1215A because the stack only includes in- 
tervening LLSs, i.e. those that are offset by at least 3000 
km s _1 from the quasar emission redshift. The peak's 
proximity to 1215 A and relatively narrow width, how- 
ever, reflect that most LLSs in our survey are located 
within Sz — 0.2 of the quasar emission redshift. 

The second strongest feature in the spectrum is the 
Lyman limit at the expected wavelength of 912A. Short- 
ward of the Lyman limit, one observes a non-zero flux 
that extends down to ~ 770A. We estimate the aver- 
age optical depth of the Lyman limit absorption in this 
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Fig. 16. — The average absorption spectrum of T912 > 2 LLSs at z as 3.7. Absorption from the H I Lyman series is readily apparent and 
one also notes strong Lyman limit absorption, as marked by the dashed vertical line. The measured average optical depth at the Lyman 
limit is T912 ~ 3.5, a value that is inconsistent (too low) with our derived f(Nm,X) distribution. We infer that the stack includes the 
contribution from a small, but non-negligible number of LLSs having rgi2 < 2. The figure also identifies a series of metal-line transitions 
arising from low and high-ions. The presence of strong O VI absorption is especially notable and suggests that LLSs are comprised of gas 
in multiple phases. 



stacked spectrum by assuming the flux at 980A provides 
a rough estimate of the absorbed continuum at the Ly- 
man limit and measure 



T912 



= - In 



/(900A) 
/(980A) 



= 3.5 



(24) 



We estimate a 20% error in this value due to effects re- 
lated to sky subtraction, uncertainty in the absorbed con- 
tinuum /(980A), and the flux- weighted average of our 
stack. 

The value of T912 may be compared to the average op- 
tical depth derived from our f(Nja,X) distribution: 



T912 



In 



J exp[-N mt T LL ]f(N m ,X)dN m 
J f(N m ,X)dN m 



(25) 



where the integrals are evaluated over the interval iVui = 
[10 17 - 5 , 10 22 ] cm -2 . Evaluating at z = 3.7 using the 
f(Nm,X) distribution function shown in Figure 10, we 
derive f%\% rs 5.4. The frequency distribution in the LLS 
regime is sufficiently flat that the higher A/hi systems 
(SLLS, DLAs) contribute significantly to the average. 

The offset between these two evaluations is most likely 
due to the inclusion of a non-negligible number of sys- 



tems having T912 < 2. Our analysis of mock spec- 
tra and our internal comparison of the T912 estimates 
for the LLS indicate that this occurs frequently. As 
described in § 5, we estimate a 0.2 dex uncertainty in 
the Niu values of LLSs with r 91 2 « 2. Although our 
tests also suggest this does not significantly affect the 
estimate of £ T >2(z), it can have a significant effect on 
the fgi2 value in the stacked spectrum. We have re- 
peated our calculation of fgi2 extending the lower limit 
of the A/hi distribution to Nm — 10 172 cm~ 2 instead of 
A/hi = 10 17 - 5 cm" 2 . For our favored f(N m ,X) distri- 
bution, we calculate T912 =4.1 and when allowing for 
sample variance with a bootstrap analysis we find con- 
sisteny in a non-negliglbe fraction of the trials (> 5%). 

Returning to the stacked spectrum (Figure 16), we 
note a series of absorption lines corresponding to strong 
metal-line transitions of low and high ionization ions. 
This suggests a highly ionized, and possibly multi-phase 
gas. The strong absorption of O VI, in particular, sug- 
gests a multi-phase medium consisting of at least one 
"cool" (T w 10 4 K), photoionized phase and another, 
more highly ionized phase which is presumably "warmer" 
(T > 10 5 K) and possibly collisionally ionized. This 
highly ionized phase has been detected in the damped 
Lya systems (Wolfe & Prochaska 2000; Fox et al. 2007a) 
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and SLLSs (Fox et al. 2007b) and its presence in our av- 
erage spectrum suggests it likely exists in lower JVhi sys- 
tems too (e.g. Simcoe et al. 2004). A more quantitative 
analysis of the metal-line absorption in LLSs, however, 
awaits high-resolution spectroscopy (e.g. Prochter et al. 
2009). 

9.4. The Physical Nature of the LLSs 

As described in the introduction, observations and 
numerical simulations associate the DLAs with high z 
galaxies residing in virializcd dark matter halos (e.g. 
M0ller et al. 2002; Pontzcn et al. 2008). The major- 
ity of absorption lines comprising the Lya forest, mean- 
while, are believed to trace Mpc-scale overdensities in 
the medium between such galaxies (e.g. Miralda-Escude 
et al. 1996). By inference, one may associate LLSs with 
lower A^hi (< 10 20 cm~ 2 ) with the interface between the 
IGM and galaxies. This inference, however, has not yet 
been extensively tested by cosmological simulations or 
empirical observation. Early works on the topic generally 
yielded too few LLSs in cosmological volumes (Katz et al. 
1996; Gardner et al. 2001; Mailer et al. 2003). More re- 
cently, Kohler & Gnedin (2007) examined LLSs in a suite 
of simulations tuned to match the observed incidence of 
LLSs at z ~ 4. Their simulations suggest LLSs are highly 
ionized gas occupying volumes of space with dimension 
1 — 100 kpc and physically associated with galaxies of a 
wide range in mass. The simulations, however, were not 
rigorously tested against observations nor did they have 
sufficient spatial resolution to "establish the physical na- 
ture of these systems" . The question remains: what is 
the physical nature of the LLSs? 

Our observations place new constraints on the struc- 
tures that give rise to LLSs absorption. The most in- 
formative measurements are £ T > 2 (X) and the shape of 
f(Nm,X) in the LLSs regime. Consider first £ T > 2 (X), 
which is proportional to the comoving number density 
?ills of the structures times their average physical size 
olls- In Figure 9, we present the t(X) values (in cu- 
mulative form) for DLAs (Prochaska & Wolfe 2009) and 
SLLSs (O'Meara et al. 2007). For the latter, we assume 
^slls(A) = 0.20 at all redshifts. We have adopted a 20% 
lower (la) value than reported by O'Meara et al. (2007) 
to crudely correct for the SDSS quasar targeting bias 
that will affect their measurement (PWO09). The figure 
demonstrates that the DLAs (especially) and the SLLSs 
have modest contributions to £ T > 2 (X). At z = 3.4, 
they contribute roughly half of the observed incidence 
of T912 > 2 LLS decreasing to w 33% by z = 4. This lat- 
ter conclusion hinges on our assumptions for £slls(X), 
in particular at z w 4 where the value is not well con- 
strained, but we expect the SLLSs to behave similarly to 
the DLAs whose incidence is not increasing significantly 
at these redshifts. We conclude that the incidence of 
T 9i2 S 10 LLSs £lls( 2 ) is comparable to that of the 
DLAs and SLLSs. As a result, it is reasonable to asso- 
ciate all LLSs with the same structures, i.e. gas within 
virialized halos as suggested by Kohler & Gnedin (2007). 
The principle challenge to this association is whether dif- 
fuse halo gas has sufficiently high cross-section to LLS ab- 
sorption. In particular, one should consider whether the 
"cold-flows" identified in numerical simulations of galaxy 
formation (Keres et al. 2005; Dekel et al. 2009) have suf- 
ficient density and size to explain the majority of LLSs. 



We are currently pursuing such analysis. 

Now consider the evolution in £ T > 2 (X) with redshift. 
As noted in § 8.2, l T > 2 (X) is observed to decrease with 
decreasing redshift. Examining Figure 9 it is evident 
that this decrease is driven by LLS with lower T912, i.e. 
by a significant decrease in £lls(A). One possible expla- 
nation for this decrease would be an increase in the in- 
tensity of the EUVB with decreasing redshift, most likely 
due to higher cmissivity from the quasar population. Un- 
like the DLAs and SLLSs systems which have T912 > 100 
and therefore have regions that are self-shielded from the 
EUVB, the T912 < 10 LLSs are probably highly ionized 
throughout (e.g. Prochter et al. 2009). Similar to the 
Lya forest, the T912 < 10 LLSs are sensitive to changes 
in the EUVB. This interpretation would predict that the 
average U parameter for T912 < 10 LLSs increase with de- 
creasing redshift. Another possible interpretation arises 
from linking LLS absorption to the cold-flows in galactic 
halos predicted by numerical simulations. The simula- 
tions reveal that the incidence of cold-flows in massive 
halos M > 1O 12 M decreases with redshift Dekel et al. 
(2009), which could result in a lower l T > 2 (X) for the 
LLSs. These assertions warrant further study and are 
likely to impact our understanding of the evolution in 
the mean free path of the universe (PWO09). 

Finally, consider the constraints on ./lls^hi, A) set 
by the observations. As described in Figure 10, we find 
that f(Nui,X) flattens at Ahi < 10 19 cm~ 2 to a power- 
law shallower than (3 « — 1. The direct implication is 
that the cross-section of gas with Ahi = 10 19 cm -2 ex- 
ceeds that of gas with Nm = 10 17 5 cm -2 . This is a 
remarkable result. Collapsed structures are generally ob- 
served to have a density gradient where the highest den- 
sity regions occupy a smaller projected cross-section than 
lower density regions. Our observations indicate the op- 
posite is true for LLSs with T912 < 10. This suggests that 
LLSs arise in structures with an extended, higher surface 
density region surrounded by a thin 'layer' of gas with 
lower Nm. Such a description brings to mind the cold- 
flows of accreting gas found in cosmological simulations 
of high z galaxies (Keres et al. 2005; Dekel et al. 2009). 
Again, we plan to explore whether such gas can explain 
the observed normalization and shape of f(Nm,X) at 
z > 3.5. 

10. SUMMARY 

In this paper, we have performed a survey for T912 > 2 
LLS absorption in the quasar spectra of the Sloan Digi- 
tal Sky Survey, Data Release 7. We established a spec- 
tral sample for statistical analysis (§ 3), estimated an 
absorbed continuum for each quasar (§ 4), searched for 
LLSs using automated algorithms (§ 5), defined the sur- 
vey path with strict criteria (§ 6), and explored the effects 
of systematic bias and uncertainty with mock spectra 
(§ 7). The primary results of this work are as follows: 

1. We measure the incidence of T912 > 2 LLSs £ T > 2 (z) 
at z = 3.3 — 4.4 and find it is well modeled as a sin- 
gle power-law, £ T > 2 (z) =C LLS [(l + z)/(l + z*)] 7LLS , 
with z* = 3.7, Clls = 1-9 ± 0.2, and 7 LLS = 
5.2 ±1.5 (68% c.L). 

2. A survey of LLSs in the SDSS spectra z cm < 
3.6 quasars confirms a previously identified bias 
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(PWO09) in the SDSS quasar targeting criteria 
that biases the sample toward sightlines with fore- 
ground LLS absorption. 

3. The number of T912 > 2 LLS per unit absorption 
length l T >2{X) is observed to decrease by w 50% 
from z = 4 to 3.4. This indicates a decrease in the 
number of systems per comoving Mpc 3 and/or a 
decrease in the average physical cross-section per 
system. We suggest it is the latter effect, possibly 
related to an increase in the EUVB with decreasing 
rcdshift or a rising radiation field local to LLSs. 

4. The measured £ T >2{X) values place an integral 
constraint on the H I frequency distribution 
f(Nm,X) at z k 3.7. Adopting previous esti- 
mates of /(Afci, X) for N m > 10 19 cm" 2 (O'Meara 
et al. 2007; Prochaska & Wolfe 2009), we constrain 
f(N m ,X) for N m = 10 17 - 5 - 10 19 cm- 2 assuming 

a power- law form f hIjS (N m ,X) = fc LLS A^ LS to 
have fc LLS w 10~ 4 - 5 and /3 LLS = -0.8 ± 0.3. This 
indicates a further shallowing of the slope as one 
decreases Nm below 10 19 cm~ 2 . 

5. Adopting constraints from the mean free path 
(PWO09) and Lya forest, we derived new 
constraints on f(Nm,X) at z « 3.7 for 
N m « 10 15 - 10 18 cm- 2 . We find that (3 = 
dlnf(Nm,X)/dlnNm must be steeper than j3 — 
-1.5 at N m w 10 15 cm- 2 . 

6. We surveyed the spectra for proximate LLSs 
(PLLSs), those with rcdshifts that arc within 
3000km s _1 of the quasar. We measure an w 25% 
lower incidence of PLLSs than intervening systems 
at z > 3.5. This lends further support to the asser- 
tion that quasars have anisotropic emission (Hcn- 
nawi & Prochaska 2007). 

Compared to previous work, our estimates of £(z) show 
systematically lower values. We suggest that the differ- 
ence is due to sample variance and/or unidentified sys- 
tematic bias in the prior analysis. This conclusion is sup- 
ported by measurements of the mean free path (PWO09) 
which do not allow for a significantly higher incidence of 
LLSs. We also find that the range of power-laws that de- 
scribe our results at z « 3.5 — 4 do not extrapolate to the 
results from z < 1 observations (Stengler-Larrea et al. 
1995). We infer that the incidence of LLSs exhibits a 
break at z « 2, qualitatively similar to that observed for 
the Lya forest (e.g. Weymann et al. 1998). The declin- 
ing incidence of LLSs per absorption length and the very 
shallow slope of f(Nm, X) at Nm < 10 19 cm -2 suggest 
that T912 < 10 LLSs arise in flattened (e.g. filamentary) 
structures that have relatively sharp edges. We associate 
these structures to the virialized halos that presumably 
give rise to SLLS and DLA absorption. Finally, we en- 
courage future work on whether such structures are con- 



sistent with the "cold-flows" identified in numerical sim- 
ulations. 

Through detailed analysis of biases and careful sample 
selection from a large and homogeneous dataset, this pa- 
per provides the first robust estimate of the incidence of 
LLSs at high rcdshift. We note that the systematic errors 
described in § 7 likely limit the precision of any future 
£(z) estimates to the order of 20 — 30%. Nevertheless, this 
is sufficient to further explore the true evolution in £(z) 
with redshift. Programs with the Hubble Space Tele- 
scope for absorption at z < 2 and with ground-based ob- 
servatories for z > 4 are currently ongoing. Altogether, 
these projects will describe the evolution of the UV back- 
ground, the growth of structure on galactic (and larger) 
scales, and the chemical enrichment history of the uni- 
verse. 
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APPENDIX 

COMPARISONS WITH KECK+LRIS SPECTRA 

To assess uncertainties (statistical and systematic) of surveying LLSs in the SDSS quasar spectra, we obtained 
independent, higher quality spectra using the LRIS spectrometer (Oke et al. 1995) on the Keck I telescope. LRIS 
employs a dichroic to split the data into two spectral channels, each with its own camera. For our observations, we 
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employed the d560 dichroic which splits the light at « 5600 A. For the blue channel, we used the 640/4000 grism which 
provides a dispersion of 0.63 A per unbinned pixel and has a nominal wavelength coverage of 3100A < A < 5600 A. 
The blue channel data was binned by 2 in both the spatial and spectral dimensions. For the red channel, we used 
the 600/7500 grating which provides a dispersion of 1.28 A per unbinned pixel, and which was tilted to provide a 
wavelength coverage of 5600A < A < 8200A. The red channel data was unbinned. All observations were obtained 
using a 1 arc-second slit which provides an ps 4 pixel FWHM corresponding to w 290km s^ 1 and w 220km s _1 for the 
blue and red data respectively. The data were obtained in good sky conditions during a 4 night run in October 2008 
and had exposure times ranging from 300 to 500 seconds. The data were reduced using the LowRedux pipeline 12 which 
bias subtracts, flat fields, optimally extracts, wavelength and flux calibrates the data to produce a final ID spectrum. 

The SDSS targets for LRIS observations were chosen to sample a range of LLSs, e.g., LLSs with Tgi 2 > 2, pLLSs 
candidates, PLLSs and spectra without apparent LLSs. Furthermore, an emphasis was placed on quasars with lower 
S/N SDSS spectra to assess the completeness of recovering LLSs. In all cases, the LRIS spectra have sufficient S/N 
to unambiguously detect the presence of absorbers with T gl2 > 1 over the full SDSS wavelength range (i.e. A > 3800 
A) for the intervening LLS survey. Particular emphasis was given to determine what absorbed continuum S/N cutoff 
should be applied to the SDSS sample. To this end, two of the authors (JXP and JMO) independently modeled LLS 
absorption in the LRIS data and compared the results to similar analysis of the the SDSS spectra. The same codes 
were used to model LLS absorption. In Table Al, we present the results of these comparisons. In Figure Al we show 
a representative sample of the LRIS data alongside their SDSS counterparts. For the SDSS data in Figure Al, we also 
show the continuum level assigned to each spectrum. 

The LRIS/SDSS comparison illustrates that with a choice of absorbed continuum S/N > 1 we recover nearly 100% 
of the LLSs with T912 > 2 in the SDSS. We see this explicitly in Table Al, where we give the values for z s tart in 
the SDSS search, where z s tart is the redshift at which the absorbed continuum S/N crosses the value of 1. In some 
cases, we identify LLSs at redshifts lower than z s t a rt m both the LRIS and SDSS data. Although these LLSs will not 
contribute to our results, they lend additional confidence in our absorbed continuum S/N cutoff. The only exceptions 
to the identification of LLSs in the LRIS and SDSS spectra at S/N > 1 are for systems with T912 ~ 2. At this optical 
depth, JXP and JMO did not always agree, even in the LRIS results. This highlights the fact that our results have an 
inherent uncertainty in log Am of w 0.2 dex. Most importantly, this disagreement docs not appear to depend on the 
S/N of the spectrum for the range of S/N we could expect from the SDSS data, and thus does not effect our choice 
for the S/N threshold. 
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SDSS-DR7 QUASARS WITH KECK+LRIS COMPARISON SPECTRA 
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Fig. Al. — Comparisons of SDSS (left panel in each pair) and follow-up Keck/LRIS observations (right panel in each pair). Overplotted 
on the SDSS spectrum is our model of the absorbed continuum (blue) and the modeled foreground Lyman limit absorption (green). For 
the SDSS spectra with S/N^ 12 > 2, we find excellent agreement for analyses of each pair of spectra. 



